Northeast Trawl Advisory Panel 2022 Restrictor Rope Experiment
Loading data setting options
Here I’m connecting to the MS Access database file. The path is wonky, but you can set it to where your file lives. I also look at the tables that are in there right now. Not 100% sure what each captures, but I tried some explorations below.
knitr::opts_chunk$set(fig.path='figs/',echo=TRUE, warning=FALSE, messages=FALSE, results = TRUE, dpi=500)
# Set specs for knitting to a github_document
## Import RODBC package (and others)
library(RODBC)
library(tidyverse)
library(lubridate)
## Set up driver info and database path
DRIVERINFO <- "Driver={Microsoft Access Driver (*.mdb, *.accdb)};"
#THIS PATH WILL NEED TO BE CHANGED
MDBPATH <- "data/NEAMAP_Restrictor_Exp.mdb"
PATH <- paste0(DRIVERINFO, "DBQ=", MDBPATH)
## Establish connection
channel <- odbcDriverConnect(PATH)
##Look at tables
tbls <- sqlTables(channel)
tbls %>% filter(TABLE_TYPE=='TABLE') %>% dplyr::select(TABLE_NAME)## TABLE_NAME
## 1 Catch
## 2 CatchNew
## 3 ChloroNew
## 4 Chlorophyll
## 5 Comments
## 6 CommentsNew
## 7 Crew
## 8 CrewNew
## 9 Cruise
## 10 CruiseNew
## 11 DogfishpupLensNew
## 12 DogfishPupsNew
## 13 Habitat
## 14 HabitatNew
## 15 Individual
## 16 IndividualNew
## 17 Inverts
## 18 LobsterNew
## 19 LongfinNew
## 20 LW_Check_Individuals
## 21 LW_Check_PanAverages
## 22 Pan
## 23 PanNew
## 24 PotTrapped
## 25 PottrappedNew
## 26 SpCodes
## 27 Station
## 28 StationNew
## 29 TagsNew
## 30 TrackData
## 31 TrackDataNew
## 32 Weather
## 33 WeatherNew
## 34 WQ
## 35 WQNew## Load data into R dataframe
Catch <- sqlQuery(channel,
"SELECT *
FROM catch;",
stringsAsFactors = FALSE)
## Load data into R dataframe (different function)
StationNew <- sqlFetch(channel, 'StationNew')
Station <- sqlFetch(channel, 'Station')
SpCodes <- sqlFetch(channel, 'SpCodes')
Comments <- sqlFetch(channel, 'Comments')
CommentsNew <- sqlFetch(channel, 'CommentsNew')
Individual <- sqlFetch(channel, 'Individual')
Pan <- sqlFetch(channel, 'Pan')
PanNew <- sqlFetch(channel, 'PanNew')
Weather <- sqlFetch(channel, 'Weather')
WQ <- sqlFetch(channel, 'WQ')
## Close and remove channel
# close(channel)
# rm(channel)
#Pulling in the species table from excel
#THIS PATH WILL NEED TO BE CHANGED
library(readxl)
NM_Species_Code_List_Color_Coordinated <- read_excel("data/NM Species Code List Color Coordinated.xls") %>% mutate(SPCODE=as.character(as.numeric(SPCODE)))A species subset
This is one way to subset the species. Summarizing the largest counts and weights and then making a vector of the species on both lists. this is not used later on, but was useful in exploring the data initially.
#making a vector of common species
common_spp_wt <- Catch %>% dplyr::group_by(VIMSCODE) %>% dplyr::summarise(TotWght=sum(TotWght)) %>% dplyr::arrange(-TotWght) %>% head(15) %>% .$VIMSCODE
common_spp_ct <- Catch %>% dplyr::group_by(VIMSCODE) %>% dplyr::summarise(COUNT=sum(COUNT)) %>% dplyr::arrange(-COUNT) %>% head(15) %>% .$VIMSCODE
common_spp <- intersect(common_spp_wt,common_spp_ct)
#common_spp
#Looking at common names of these species
NM_Species_Code_List_Color_Coordinated %>% filter(SPCODE %in% common_spp) %>% print()## # A tibble: 10 × 13
## SPCODE COMMON LATIN PHYLUM CLASS ORDER FAMILY GROUP HABIT PROCESS MERISTIC
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 37 Atlantic… Brev… Chord… Acti… Clup… Clupe… bony… <NA> N <NA>
## 2 4 butterfi… Pepr… Chord… Acti… Perc… Strom… bony… <NA> <NA> <NA>
## 3 71 northern… Prio… Chord… Acti… Scor… Trigl… bony… <NA> <NA> <NA>
## 4 15 red hake Urop… Chord… Acti… Gadi… Gadid… bony… <NA> N <NA>
## 5 1 scup Sten… Chord… Acti… Perc… Spari… bony… <NA> <NA> <NA>
## 6 6 silver h… Merl… Chord… Acti… Gadi… Gadid… bony… <NA> <NA> <NA>
## 7 44 winter f… Pseu… Chord… Acti… Pleu… Pleur… bony… <NA> N <NA>
## 8 171 little s… Leuc… Chord… Elas… Raji… Rajid… cart… <NA> M TL
## 9 167 spiny do… Squa… Chord… Elas… Squa… Squal… cart… <NA> M TL
## 10 1028 longfin … Loli… Mollu… Ceph… Teut… Lolig… ceph… <NA> <NA> <NA>
## # ℹ 2 more variables: `date added` <dttm>, COMMENT <chr>Catch %>% dplyr::group_by(VIMSCODE) %>% dplyr::summarise(N_STATIONS=n_distinct(STATION)) %>% dplyr::arrange(-N_STATIONS)## # A tibble: 73 × 2
## VIMSCODE N_STATIONS
## <int> <int>
## 1 1028 140
## 2 4 135
## 3 1 134
## 4 171 131
## 5 3 120
## 6 163 119
## 7 44 118
## 8 6 110
## 9 2 85
## 10 15 83
## # ℹ 63 more rowslibrary(patchwork)
#Making a plot of the numbers weights
a <- Catch %>% dplyr::group_by(VIMSCODE) %>% dplyr::summarise(AvgTotWght=mean(TotWght),TotWght=sum(TotWght)) %>% dplyr::arrange(-TotWght) %>% left_join(.,
Catch %>% dplyr::group_by(VIMSCODE) %>% dplyr::summarise(AvgTotCOUNT=mean(COUNT),COUNT=sum(COUNT)) %>% dplyr::arrange(-COUNT)) %>%
left_join(.,Catch %>% dplyr::group_by(VIMSCODE) %>% dplyr::summarise(N_STATIONS=n_distinct(STATION)) %>% dplyr::arrange(-N_STATIONS)) %>%
left_join(.,NM_Species_Code_List_Color_Coordinated %>% mutate(VIMSCODE=as.numeric(SPCODE)) %>%
dplyr::select(VIMSCODE,COMMON),by=('VIMSCODE')) %>%
mutate(PER_STAT = round(1-N_STATIONS/max(N_STATIONS+1),2)) %>% arrange(PER_STAT) %>%
head(35) %>%
pivot_longer(names_to = 'metric',values_to = 'value',-COMMON) %>%
filter(metric != 'VIMSCODE') %>%
mutate(COMMON = fct_reorder(COMMON, value)) %>%
ggplot() + geom_bar(aes(x=COMMON,y=value),stat='identity') + facet_wrap(~metric,scales='free') + coord_flip() + scale_y_log10()
b <- Catch %>% dplyr::group_by(VIMSCODE) %>% dplyr::summarise(TotWght=sum(TotWght)) %>% dplyr::arrange(-TotWght) %>% left_join(.,
Catch %>% dplyr::group_by(VIMSCODE) %>% dplyr::summarise(COUNT=sum(COUNT)) %>% dplyr::arrange(-COUNT)) %>%
left_join(.,Catch %>% dplyr::group_by(VIMSCODE) %>% dplyr::summarise(N_STATIONS=n_distinct(STATION)) %>% dplyr::arrange(-N_STATIONS)) %>%
left_join(.,NM_Species_Code_List_Color_Coordinated %>% mutate(VIMSCODE=as.numeric(SPCODE)) %>%
dplyr::select(VIMSCODE,COMMON),by=('VIMSCODE')) %>%
pivot_longer(names_to = 'metric',values_to = 'value',-COMMON) %>%
filter(metric != 'VIMSCODE') %>%
mutate(COMMON = fct_reorder(COMMON, value)) %>%
head(35) %>%
filter(metric == 'COUNT',!COMMON %in% c('butterfish','scup','longfin inshore squid','Atlantic menhaden','silver hake')) %>%
ggplot() + geom_bar(aes(x=COMMON,y=value),stat='identity') + facet_wrap(~metric,scales='free') + coord_flip()
c <- Catch %>%
left_join(.,Catch %>% dplyr::group_by(VIMSCODE) %>% dplyr::summarise(N_STATIONS=n_distinct(STATION)) %>% dplyr::arrange(-N_STATIONS)) %>%
left_join(.,NM_Species_Code_List_Color_Coordinated %>% mutate(VIMSCODE=as.numeric(SPCODE)) %>%
dplyr::select(VIMSCODE,COMMON),by=('VIMSCODE')) %>%
#head(35) %>%
filter(!VIMSCODE %in% c(1,4,1028)) %>%
ggplot() + geom_boxplot(data= . %>% filter(COMMON != 'Other'),aes(x=fct_lump(fct_reorder(COMMON,COUNT),15),y=COUNT)) + coord_flip() + scale_y_log10()
d <- Catch %>%
left_join(.,Catch %>% dplyr::group_by(VIMSCODE) %>% dplyr::summarise(N_STATIONS=n_distinct(STATION)) %>% dplyr::arrange(-N_STATIONS)) %>%
left_join(.,NM_Species_Code_List_Color_Coordinated %>% mutate(VIMSCODE=as.numeric(SPCODE)) %>%
dplyr::select(VIMSCODE,COMMON),by=('VIMSCODE')) %>%
#head(35) %>%
filter(!VIMSCODE %in% c(1,4,1028)) %>%
ggplot() + geom_boxplot(data= . %>% filter(COMMON != 'Other'),aes(x=fct_lump(fct_reorder(COMMON,TotWght),15),y=TotWght)) + coord_flip() + scale_y_log10()Ordination of stations
I’m using the lat 4 digits of the station code to plot out the catches in space. The species values are shown in red to highlight the specific catch elements that are driving patterns. I was thinking about this as maybe a way to double check for stations or pairs of stations that are outliers (0424 and 9424 for example). Right now it’s just on the ~20 species we identify above.
library(vegan)
library(ggvegan)
library(ggrepel)
#Looking at catch data weight
catch_ord <- Catch %>% as.tibble() %>%
filter(VIMSCODE %in% common_spp_wt) %>% #COMMENT OUT TO DO ON FULL SET OF SPECIES
inner_join(.,Station) %>% mutate(STATION=str_sub(STATION,-4)) %>%
mutate(SEASON=case_when(month(DATE) < 7 ~ 'S',TRUE~'F')) %>%
mutate(STATION=paste(SEASON,STATION,sep='')) %>%
group_by(VIMSCODE,Restrictor,STATION) %>% summarise(TotWght=sqrt(sum(TotWght))) %>%
pivot_wider(names_from = c(VIMSCODE),values_from = TotWght,values_fill = 0) %>% as_tibble()
catch_mds <- metaMDS(catch_ord %>% dplyr::select(-Restrictor,-STATION), distance = "bray", autotransform = TRUE,tidy=TRUE)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1627368
## Run 1 stress 0.163298
## Run 2 stress 0.1627369
## ... Procrustes: rmse 0.0009123008 max resid 0.009364886
## ... Similar to previous best
## Run 3 stress 0.1642373
## Run 4 stress 0.1627375
## ... Procrustes: rmse 0.001044289 max resid 0.01076081
## Run 5 stress 0.1632985
## Run 6 stress 0.164782
## Run 7 stress 0.1627391
## ... Procrustes: rmse 0.00128301 max resid 0.01328321
## Run 8 stress 0.1647814
## Run 9 stress 0.1632987
## Run 10 stress 0.1627378
## ... Procrustes: rmse 0.0002718825 max resid 0.002325157
## ... Similar to previous best
## Run 11 stress 0.1633081
## Run 12 stress 0.1642367
## Run 13 stress 0.165025
## Run 14 stress 0.1642356
## Run 15 stress 0.1643469
## Run 16 stress 0.1627373
## ... Procrustes: rmse 0.000122544 max resid 0.00126336
## ... Similar to previous best
## Run 17 stress 0.1642391
## Run 18 stress 0.1632983
## Run 19 stress 0.1642368
## Run 20 stress 0.1632983
## *** Best solution repeated 3 times# test_set <- fortify(catch_mds) %>% filter(Score == 'sites') %>% cbind(.,catch_ord$STATION) %>% dplyr::rename(STATION=`catch_ord$STATION`) %>%
# mutate(Restrictor=str_sub(STATION,1,1)) %>%
# mutate(PAIR=str_sub(STATION,-3)) %>%
# filter(abs(NMDS1) < 1.5 & abs(NMDS2) < 1) %>%
# count(PAIR) %>%
# filter(n == 2) %>% .$PAIR
#making some plots
library(patchwork)
library(ggthemes)
a <- fortify(catch_mds) %>% filter(Score == 'sites') %>%
cbind(.,catch_ord$STATION) %>% dplyr::rename(STATION=`catch_ord$STATION`) %>%
mutate(Restrictor=str_sub(STATION,2,2)) %>%
mutate(PAIR=str_sub(STATION,-3)) %>%
#filter(PAIR %in% test_set) %>%
ggplot(aes(x=NMDS1,y=NMDS2)) + geom_point() +
geom_label_repel(aes(label=STATION)) +
theme(legend.position = 'None')
b <- fortify(catch_mds) %>% filter(Score != 'sites') %>%
left_join(.,NM_Species_Code_List_Color_Coordinated %>%
#mutate(SPCODE=as.character(as.numeric(SPCODE))) %>%
dplyr::select(VIMSCODE=SPCODE,COMMON), by = c('Label'='VIMSCODE')) %>%
ggplot() +
geom_segment(aes(x=0,xend=NMDS1,y=0,yend=NMDS2),colour='Red',alpha=0.5) +
geom_label(aes(x=NMDS1,y=NMDS2,label=COMMON),colour='Red')
a + b
#Season
s <- fortify(catch_mds) %>% filter(Score == 'sites') %>%
cbind(.,catch_ord$STATION) %>% dplyr::rename(STATION=`catch_ord$STATION`) %>%
mutate(Season=str_sub(STATION,1,1)) %>%
mutate(PAIR=str_sub(STATION,-3)) %>%
#filter(PAIR %in% test_set) %>%
mutate(Season = case_when(Season == 'F' ~ 'Fall',TRUE ~ 'Spring')) %>%
mutate(Season = factor(Season,levels = c("Spring", "Fall"))) %>%
ggplot(aes(x=NMDS1,y=NMDS2,colour=Season)) + geom_point() +
#geom_label_repel(aes(label=STATION)) +
#theme(legend.position = 'None') +
theme(legend.position='top',
legend.justification='left',
legend.direction='horizontal') +
scale_color_viridis_d()
#Restrictor
r <- fortify(catch_mds) %>% filter(Score == 'sites') %>%
cbind(.,catch_ord$STATION) %>% dplyr::rename(STATION=`catch_ord$STATION`) %>%
mutate(Restrictor=str_sub(STATION,2,2)) %>%
mutate(PAIR=str_sub(STATION,-3)) %>%
#filter(PAIR %in% test_set) %>%
mutate(Restrictor = case_when(Restrictor == '9' ~ 'Yes',TRUE ~ 'No')) %>%
ggplot(aes(x=NMDS1,y=NMDS2,colour=Restrictor)) + geom_point() +
#geom_label_repel(aes(label=STATION)) +
scale_color_tableau() +
theme(legend.position='top',
legend.justification='left',
legend.direction='horizontal')
catch_ordination <- s + r
catch_ordination
Making a map of the sampling sites
#Making a quick map of stations
#Library loading
library(tidyverse)
library(marmap)## Registered S3 methods overwritten by 'adehabitatMA':
## method from
## print.SpatialPixelsDataFrame sp
## print.SpatialPixels sp##
## Attaching package: 'marmap'## The following object is masked from 'package:grDevices':
##
## as.raster#for loading the raster library
library(pals)
library(ggthemes)
library(ggnewscale)
library(sf)## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUElibrary(ggspatial)
#Read in shapefiles
US.areas <- rgdal::readOGR("./shapefiles","USA")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\andrew.jones\Desktop\PEMAD_CRB_NTAP_RESTRICTOR_ROPE_EXP\shapefiles", layer: "USA"
## with 24 features
## It has 4 fieldscanada.areas <-rgdal::readOGR("./shapefiles","Canada")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\andrew.jones\Desktop\PEMAD_CRB_NTAP_RESTRICTOR_ROPE_EXP\shapefiles", layer: "Canada"
## with 4998 features
## It has 11 fieldsatl <- marmap::getNOAA.bathy(-72,-70, 41, 42,resolution=1)## Querying NOAA database ...## This may take seconds to minutes, depending on grid size## Building bathy matrix ...#atl <- read.bathy("./bathydata/marmap_coord_-80_33_-66_46_res_1.csv",header = TRUE)
atl = fortify.bathy(atl)
library(inlabru)## Loading required package: sp#Setting the map boundary
lons = c(min(Station$LONGITUDE_START)-0.1, max(Station$LONGITUDE_START)+0.1)
lats = c(min(Station$LATITUDE_START)-0.1, max(Station$LATITUDE_START)+0.1)
#Making a function to bin the catches
label_interval <- function(breaks) {
paste0("(", breaks[1:length(breaks) - 1], "-", breaks[2:length(breaks)], ")")
}
# Select for start geometries and convert to sf object
pnts_start <- Station %>% st_as_sf( coords = c("LONGITUDE_START", "LATITUDE_START"), crs = 4326)
# Select for end geometries and convert to sf object
pnts_end <- Station %>% st_as_sf(coords = c("LONGITUDE_END", "LATITUDE_END"), crs = 4326)
# Combine start and end geometries
cbind(pnts_start,pnts_end) -> points_ready
# Generate line segments via union of paired geometries
line_segments <- as.data.frame(st_sfc(mapply(function(a,b){st_cast(st_union(a,b),"LINESTRING")}, points_ready$geometry, points_ready$geometry.1, SIMPLIFY=FALSE)))
line_segments <- cbind(line_segments,Station) %>% mutate(PAIR=str_sub(STATION,-3)) %>%
mutate(SEASON=case_when(month(DATE) < 7 ~ 'Spring',TRUE~'Fall')) %>%
mutate(PAIR=paste(SEASON,PAIR,sep=''))
line_segments$SEASON <- factor(line_segments$SEASON,
levels = c("Spring", "Fall"))
line_segments <- line_segments %>% st_as_sf(crs = 4326)
my_depths <- c(0,10,20,30,40,50,60)
sampling_map <- ggplot() +
#First bathymetry
geom_contour_filled(data=atl,
aes(x=x,y=y,z=-1*z),
breaks=c(-50,10,20,30,40,50,60),
size=c(0.1),alpha=0.75) +
geom_sf() +
scale_fill_manual(values = brewer.blues(10)[5:10],name=paste("Depth (m)"),labels = label_interval(my_depths)) +
theme(legend.position = 'bottom', legend.box="vertical", legend.margin=margin()) +
guides(fill=guide_legend(nrow=2)) +
theme(legend.key.width=unit(0.25, "cm")) +
new_scale_fill() +
#Now the catch data
#Canada
geom_sf(data=canada.areas %>% st_as_sf() %>% st_transform(., 4326)) +
#The US
geom_sf(data=US.areas %>% st_as_sf() %>% st_transform(., 4326)) +
geom_sf(data = line_segments, aes(colour=SEASON), size=1) +
coord_sf(xlim = lons, ylim = lats, crs="+proj=longlat +datum=WGS84") +
theme_map() +
#theme(legend.position = 'bottom', legend.box="horizontal", legend.margin=margin()) +
theme(legend.position = c(1, 0), legend.justification = c(1, 0),
legend.spacing.y = unit(0.2, "mm"),
#panel.border = element_rect(colour = "black", fill=NA),
aspect.ratio = 1, axis.text = element_text(colour = 1, size = 14),
legend.background = element_rect(colour = 'White',fill='White'),
legend.box.background = element_rect(colour = 'White',fill='White')
) +
theme(legend.key.width=unit(0.25, "cm"),legend.text=element_text(size=8)) +
scale_colour_viridis_d() +
#scale_color_manual(values = 'Midnight blue') +
guides(fill=guide_legend(nrow=2),color=guide_legend(nrow=1),size=guide_legend(nrow=2)) +
annotation_scale(location = "tl", width_hint = 0.3) +
annotation_north_arrow(location = "tl", which_north = "true",
pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
style = north_arrow_fancy_orienteering) +
labs(colour='Season')Adding the solar zenith info
library(fishmethods)
zen_data <- Station %>% as_tibble() %>% mutate(hour=hour(TowBegin),day=day(TowBegin),month=month(TowBegin),year=year(TowBegin),timezone = -5)
station_zen_data <- zen_data %>% cbind(.,astrocalc4r(hour=zen_data$hour,day=zen_data$day,month=zen_data$month,
year=zen_data$year,timezone = zen_data$timezone,
lat = zen_data$LATITUDE_END,lon = zen_data$LONGITUDE_END,
withinput = FALSE, seaorland = "maritime", acknowledgment = FALSE))
station_zen_data %>% dplyr::select(month,zenith,hour) %>% ggplot(aes(x=hour,y=zenith,colour=month)) + geom_point()
pair_zen_data <- station_zen_data %>%
mutate(PAIR=str_sub(STATION,-3)) %>%
mutate(SEASON=case_when(month(DATE) < 7 ~ 'S',TRUE~'F')) %>%
mutate(PAIR=paste(SEASON,PAIR,sep='')) %>%
group_by(PAIR,DATE) %>% summarise(zenith=mean(zenith)) %>%
mutate(DATE=as.POSIXct(DATE,format = '%Y-%M-%D')) %>%
mutate(SEASON=case_when(month(DATE) < 7 ~ 'Spring',TRUE~'Fall'))## `summarise()` has grouped output by 'PAIR'. You can override using the
## `.groups` argument.Looking at the distance between tows
#Looking at the distance between tows
restric_y <- line_segments %>% select(STATION,DATE,geometry,Restrictor) %>% mutate(PAIR=str_sub(STATION,-3)) %>%
mutate(SEASON=case_when(month(DATE) < 7 ~ 'S',TRUE~'F')) %>% mutate(PAIR=paste(SEASON,PAIR,sep='')) %>%
dplyr::select(PAIR,Restrictor,geometry) %>% filter(Restrictor == 'Y') %>% rename(geometry_y = geometry)
restric_n <- line_segments %>% select(STATION,DATE,geometry,Restrictor) %>% mutate(PAIR=str_sub(STATION,-3)) %>%
mutate(SEASON=case_when(month(DATE) < 7 ~ 'S',TRUE~'F')) %>% mutate(PAIR=paste(SEASON,PAIR,sep='')) %>%
dplyr::select(PAIR,Restrictor,geometry) %>% filter(Restrictor != 'Y') %>% rename(geometry_n = geometry)
#ggplot(data=restric_n,aes(x=as.numeric(str_sub(restric_y$PAIR,2,4)),y=as.numeric(str_sub(restric_n$PAIR,2,4)))) + geom_point()
#Looking at the tows
num <- 8
ggplot() +
geom_sf(data=restric_y$geometry_y) +
geom_sf(data=restric_n$geometry_n,colour='red')
#1/4 mile is 400 meters
st_distance(restric_y$geometry_y[num],restric_n$geometry_n[num])## Units: [m]
## [,1]
## [1,] 205.532723#Mean distance
dist <- st_distance(restric_y$geometry_y,restric_n$geometry_n,by_element = FALSE) %>% diag() %>% cbind(.,as.character(str_sub(restric_y$PAIR,1,4))) %>% as.data.frame()
colnames(dist) <- c("Dist. (m)", "PAIR")
#Min distance (number crossing)
dist <- st_distance(restric_y$geometry_y,restric_n$geometry_n,by_element = FALSE) %>% diag() %>% as.numeric() %>% as_tibble() %>% cbind(.,as.character(str_sub(restric_y$PAIR,1,4))) %>% as.tibble()
colnames(dist) <- c("Dist. (m)", "PAIR")
dist %>%
ggplot(aes(x=`Dist. (m)`)) + geom_histogram() + geom_vline(xintercept = 400)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#Looking at which PAIR was far apart
dist %>% filter(`Dist. (m)` > 400)## # A tibble: 3 × 2
## `Dist. (m)` PAIR
## <dbl> <chr>
## 1 1562. S265
## 2 603. F149
## 3 536. F424#Looking at crossing tows
crossing <- dist %>% filter(`Dist. (m)` == 0) %>% .$PAIR
length(crossing)## [1] 24Adding in the current direction relative to the tow
As reported by the survey crew
#Looking current direction in pair
#Codes from Dustin's email
current_codes <- tribble(
~TOW_DIRECTION_CURRENT, ~CURR_DESC,
'1', 'Fair Tide',
'2', 'Head Tide',
'3', 'Perpendicular',
'4', 'Oblique-Fair',
'5', 'Oblique-Head',
'6', 'Slack')
station_pair_current <- Station %>% mutate(PAIR=str_sub(STATION,-3)) %>%
mutate(SEASON=case_when(month(DATE) < 7 ~ 'S',TRUE~'F')) %>%
mutate(PAIR=paste(SEASON,PAIR,sep='')) %>%
dplyr::select(Restrictor,PAIR,TOW_DIRECTION_CURRENT) %>%
left_join(.,current_codes) %>%
dplyr::select(-TOW_DIRECTION_CURRENT) %>%
pivot_wider(names_from = Restrictor,values_from = CURR_DESC) %>%
mutate(CURRENT=paste(N,'->',Y)) %>%
dplyr::select(PAIR,CURRENT) %>%
mutate(CURRENT_DIF = case_when(
CURRENT == 'Head Tide -> Head Tide' ~ 'SAME',
CURRENT == 'Perpendicular -> Perpendicular' ~ 'SAME',
CURRENT == 'Fair Tide -> Fair Tide' ~ 'SAME',
TRUE ~ 'DIFFERENT'
))## Joining with `by = join_by(TOW_DIRECTION_CURRENT)`#Looking at how prevalent each of these was
station_pair_current %>% group_by(CURRENT,CURRENT_DIF) %>% tally() %>% arrange(-n)## # A tibble: 12 × 3
## # Groups: CURRENT [12]
## CURRENT CURRENT_DIF n
## <chr> <chr> <int>
## 1 Head Tide -> Head Tide SAME 30
## 2 Perpendicular -> Perpendicular SAME 16
## 3 Fair Tide -> Fair Tide SAME 13
## 4 Head Tide -> Fair Tide DIFFERENT 3
## 5 Slack -> Head Tide DIFFERENT 2
## 6 Fair Tide -> Head Tide DIFFERENT 1
## 7 Fair Tide -> Perpendicular DIFFERENT 1
## 8 Head Tide -> Perpendicular DIFFERENT 1
## 9 NA -> Perpendicular DIFFERENT 1
## 10 Perpendicular -> Head Tide DIFFERENT 1
## 11 Perpendicular -> NA DIFFERENT 1
## 12 Slack -> Fair Tide DIFFERENT 1Assembling Aggregate Weight Data
Putting the catch data together. Comparing weight estimates accounting for differences in swept area and not accounting for differences in swept area. Swept area seems to have little impact on weights.
#Trying out some linear models
catch_rest <- Catch %>% as.tibble() %>%
complete(
STATION,VIMSCODE,
fill = list(TotWght = 0,COUNT=0)
) %>%
inner_join(.,Station) %>%
mutate(SA=TowDistTrack*NetWidth,MEAN_SA = mean(SA),
Adj_SA = SA/MEAN_SA,TotWght_Adj=TotWght*(Adj_SA)) %>%
mutate(PAIR=str_sub(STATION,-3)) %>%
mutate(SEASON=case_when(month(DATE) < 7 ~ 'S',TRUE~'F')) %>%
mutate(PAIR=paste(SEASON,PAIR,sep='')) %>%
group_by(VIMSCODE,Restrictor,PAIR) %>% summarise(TotWght=sum(TotWght)) %>%
ungroup() %>%
pivot_wider(names_from = c(Restrictor,VIMSCODE),
values_from = TotWght,values_fill = 0) %>%
as_tibble() %>%
pivot_longer(values_to = c('TotWght'),names_to = 'REST_SPP',-PAIR) %>%
mutate(Restrictor=str_sub(REST_SPP,1,1),VIMSCODE=str_sub(REST_SPP,3)) %>%
dplyr::select(-REST_SPP) %>%
pivot_wider(names_from = Restrictor,values_from = TotWght,values_fill=0) %>%
mutate(N=log10(N+1),Y=log10(Y+1)) %>%
left_join(.,NM_Species_Code_List_Color_Coordinated %>%
dplyr::select(VIMSCODE=SPCODE,COMMON), by = c('VIMSCODE'='VIMSCODE')) %>%
left_join(.,station_pair_current %>% dplyr::select(PAIR,CURRENT_DIF)) %>%
left_join(.,pair_zen_data %>% dplyr::select(PAIR,zenith)) %>%
left_join(.,dist)
#Total pairs w/o a couple
catch_rest$PAIR %>% n_distinct()## [1] 71#Catch with and without swept area
#Checking to see if there is a need to adjust for swept area
Catch %>% as.tibble() %>%
complete(
STATION,VIMSCODE,
fill = list(TotWght = 0,COUNT=0)
) %>%
inner_join(.,Station) %>%
mutate(SA=TowDistTrack*NetWidth,MEAN_SA = mean(SA),
Adj_SA = SA/MEAN_SA,TotWght_Adj=TotWght*(Adj_SA)) %>% #CONVERTING FOR SWEPT AREA
filter(VIMSCODE %in% c(1028,4,1,171,44,6,15)) %>%
ggplot(aes(TotWght,TotWght_Adj,color=VIMSCODE)) + geom_point() + facet_wrap(~VIMSCODE,scales='free')
#Making a plot of catch in control and treatment nets relative to a 1:1 line
#Figure S1
library(ggExtra)
p <- catch_rest %>% filter(VIMSCODE %in% c(1028,4,1,171,44,6,15)) %>%
ggplot(aes(x=N,y=Y,colour=COMMON)) + geom_point() #+ facet_wrap(~COMMON)
p <- p + geom_abline(slope = 1,intercept = 0) +
labs(color='Species',fill='Species',x='Control Catch log10(kgs + 1)',
y='Treatment Catch log10(kgs + 1)') +
scale_color_tableau() + scale_fill_tableau() +
theme(legend.position="bottom")
p + theme(legend.position="none") -> raw_weight_plot
#ggMarginal(p, type="density") -> raw_weight_plot
raw_weight_plot
ind_sp_raw_weight_plot <- p + geom_abline(slope = 1,intercept = 0) +
labs(color='Species',fill='Species',x='Control Catch log10(kgs + 1)',
y='Treatment Catch log10(kgs + 1)') +
scale_color_tableau() + scale_fill_tableau() +
theme(legend.position="none") + facet_wrap(~COMMON,scales = 'free')
raw_weight_plot + ind_sp_raw_weight_plot + plot_annotation(tag_levels = 'A')
A table of catch statisitcs
This table becomes Table 2 and helps us identify the species of interest for the study.
#Making a table of catch stats
catch_table <- Catch %>% as.tibble() %>%
complete(
STATION,VIMSCODE,
fill = list(TotWght = 0,COUNT=0)) %>%
inner_join(.,Station) %>%
mutate(PAIR=str_sub(STATION,-3)) %>%
mutate(SEASON=case_when(month(DATE) < 7 ~ 'S',TRUE~'F')) %>%
mutate(PAIR=paste(SEASON,PAIR,sep='')) %>%
mutate(POS_STATION=case_when(TotWght > 1 ~ 1, TRUE ~ 0)) %>%
group_by(VIMSCODE,Restrictor) %>%
filter(VIMSCODE %in% c(1028,4,1,171,44,6,15)) %>%
dplyr::summarise(AvgTotWght=mean(TotWght),
SETotWght=sd(TotWght)/sqrt(142),
CVTotWght=sd(TotWght)/AvgTotWght,
TotWght=sum(TotWght),
AvgTotCOUNT=mean(COUNT),
SETotCOUNT=sd(COUNT)/sqrt(142),
CVTotCOUNT=sd(COUNT)/AvgTotCOUNT,
COUNT=sum(COUNT),
N_STATIONS=sum(POS_STATION)) %>%
dplyr::arrange(-N_STATIONS) %>%
left_join(.,NM_Species_Code_List_Color_Coordinated %>% mutate(VIMSCODE=as.numeric(SPCODE)) %>%
dplyr::select(VIMSCODE,COMMON),by=('VIMSCODE')) %>%
mutate(PER_STAT = round(N_STATIONS/max(N_STATIONS+2),2)) %>%
arrange(PER_STAT) %>% arrange(desc(N_STATIONS)) %>%
filter(AvgTotCOUNT > 10,AvgTotWght>2.5,PER_STAT > 0.5) %>%
relocate('Restrictor'=Restrictor,'Species'=COMMON,
'Total positive stations'=N_STATIONS,
'Proportion positive stations'=PER_STAT,
'Total weight (kgs)'=TotWght,
'Station total count'=COUNT,
'Station mean weight (kgs)'=AvgTotWght,
'Standard error station weight (kgs)'=SETotWght,
'Coefficient of variation of station weight (kgs)'=CVTotWght,
'Station mean count'=AvgTotCOUNT,
'Standard error of station count'=SETotCOUNT,
'Coefficient of variation of station count'=CVTotCOUNT)## Joining with `by = join_by(STATION)`
## `summarise()` has grouped output by 'VIMSCODE'. You can override using the
## `.groups` argument.library(gt)
ct <- catch_table %>% ungroup() %>% dplyr::select(-VIMSCODE) %>%
gt(groupname_col = c("Species"),rowname_col = "Restrictor") %>%
fmt_number(columns=c(4),decimals = c(2)) %>%
fmt_number(columns=c(5),decimals = c(0)) %>%
fmt_number(columns=c(6:12),decimals = c(1)) %>%
cols_label(Restrictor = md("**Restrictor**"),
'Total positive stations' = md("**Tot. positive stations**"),
'Species' = md("**Species**"),
'Proportion positive stations' = md("**Prop. positive stations**"),
'Total weight (kgs)' = md("**Total weight (kgs)**"),
'Station total count' = md("**Station total count**"),
'Station mean weight (kgs)' = md("**Station mean weight (kgs)**"),
'Standard error station weight (kgs)' = md("**SE station weight (kgs)**"),
'Coefficient of variation of station weight (kgs)' = md("**CV of station weight (kgs)**"),
'Station mean count' = md("**Station mean count**"),
'Standard error of station count' = md("**SE of station count**"),
'Coefficient of variation of station count' = md("**CV of station count**"))
ct #%>% gtsave("tab_2.html") #uncomment to save as html | Tot. positive stations | Prop. positive stations | Total weight (kgs) | Station total count | Station mean weight (kgs) | SE station weight (kgs) | CV of station weight (kgs) | Station mean count | SE of station count | CV of station count | |
|---|---|---|---|---|---|---|---|---|---|---|
| longfin inshore squid | ||||||||||
| Y | 71 | 0.97 | 2,734 | 96,787.0 | 38.5 | 7.5 | 2.3 | 1,363.2 | 180.7 | 1.6 |
| N | 68 | 0.93 | 1,739 | 72,938.0 | 24.5 | 2.1 | 1.0 | 1,027.3 | 120.9 | 1.4 |
| little skate | ||||||||||
| Y | 64 | 0.97 | 832 | 1,710.0 | 11.7 | 0.9 | 0.9 | 24.1 | 1.9 | 0.9 |
| N | 63 | 0.95 | 910 | 1,849.0 | 12.8 | 1.1 | 1.1 | 26.0 | 2.4 | 1.1 |
| butterfish | ||||||||||
| N | 61 | 0.97 | 8,140 | 721,698.0 | 114.6 | 13.9 | 1.4 | 10,164.8 | 1,844.2 | 2.2 |
| Y | 60 | 0.95 | 11,988 | 1,015,687.0 | 168.8 | 29.8 | 2.1 | 14,305.5 | 3,746.2 | 3.1 |
| scup | ||||||||||
| N | 53 | 0.96 | 3,074 | 214,630.0 | 43.3 | 6.1 | 1.7 | 3,023.0 | 545.1 | 2.1 |
| Y | 48 | 0.87 | 3,178 | 259,765.0 | 44.8 | 5.6 | 1.5 | 3,658.7 | 630.2 | 2.1 |
| winter flounder | ||||||||||
| Y | 47 | 0.96 | 427 | 2,374.0 | 6.0 | 0.8 | 1.6 | 33.4 | 4.4 | 1.6 |
| N | 43 | 0.88 | 413 | 2,303.0 | 5.8 | 0.8 | 1.6 | 32.4 | 3.9 | 1.4 |
| silver hake | ||||||||||
| Y | 38 | 0.95 | 548 | 6,921.0 | 7.7 | 1.9 | 2.9 | 97.5 | 23.7 | 2.9 |
| N | 37 | 0.92 | 600 | 7,288.0 | 8.4 | 1.4 | 2.0 | 102.6 | 16.2 | 1.9 |
| red hake | ||||||||||
| Y | 28 | 0.93 | 214 | 1,609.0 | 3.0 | 0.7 | 2.6 | 22.7 | 5.8 | 3.0 |
More catch composition stats
Also making station summary tables
#Looking at the number of number of paired positives
Catch %>% as.tibble() %>% mutate(VIMSCODE=as.character(VIMSCODE)) %>%
left_join(.,NM_Species_Code_List_Color_Coordinated %>%
dplyr::select(VIMSCODE=SPCODE,COMMON), by = c('VIMSCODE'='VIMSCODE')) %>%
group_by(COMMON,VIMSCODE) %>% tally() %>%
arrange(-n) %>% mutate(prop_pos=n/142) %>% head(20)## # A tibble: 20 × 4
## # Groups: COMMON [20]
## COMMON VIMSCODE n prop_pos
## <chr> <chr> <int> <dbl>
## 1 longfin inshore squid 1028 140 0.986
## 2 butterfish 4 135 0.951
## 3 scup 1 134 0.944
## 4 little skate 171 131 0.923
## 5 summer flounder 3 120 0.845
## 6 smooth dogfish 163 119 0.838
## 7 winter flounder 44 118 0.831
## 8 silver hake 6 110 0.775
## 9 black seabass 2 85 0.599
## 10 fourspot flounder 237 83 0.585
## 11 jonah crab 600 83 0.585
## 12 red hake 15 83 0.585
## 13 winter skate 173 83 0.585
## 14 Gulf Stream flounder 247 79 0.556
## 15 northern searobin 71 75 0.528
## 16 American lobster 20 67 0.472
## 17 striped searobin 72 66 0.465
## 18 American shad 30 60 0.423
## 19 windowpane 63 60 0.423
## 20 alewife 26 52 0.366station_pair_info <- Station %>%
mutate(PAIR=str_sub(STATION,-3)) %>%
mutate(SEASON=case_when(month(DATE) < 7 ~ 'S',TRUE~'F')) %>%
mutate(PAIR=paste(SEASON,PAIR,sep='')) %>%
dplyr::select(PAIR,DATE,DepthBeg) %>%
mutate(DATE=as.POSIXct(DATE,format = '%Y-%M-%D')) %>%
mutate(SEASON=case_when(month(DATE) < 7 ~ 'Spring',TRUE~'Fall')) %>%
group_by(PAIR,DATE,SEASON) %>% summarise(DEPTH = mean(DepthBeg,na.rm=TRUE))## `summarise()` has grouped output by 'PAIR', 'DATE'. You can override using the
## `.groups` argument.station_table <- Station %>%
mutate(PAIR=str_sub(STATION,-3)) %>%
mutate(SEASON=case_when(month(DATE) < 7 ~ 'S',TRUE~'F')) %>%
mutate(PAIR=paste(SEASON,PAIR,sep='')) %>%
#dplyr::select(PAIR,DATE,DepthBeg) %>%
mutate(DATE=as.POSIXct(DATE,format = '%Y-%M-%D')) %>%
mutate(SEASON=case_when(month(DATE) < 7 ~ 'Spring',TRUE~'Fall')) %>%
mutate(PAIR=str_sub(STATION,-4)) %>%
mutate(Restrictor=str_sub(PAIR,1,1)) %>%
mutate(Restrictor = case_when(Restrictor == '9' ~ 'Yes',TRUE ~ 'No')) %>%
group_by(Restrictor,SEASON) %>%
summarise(Stations = n_distinct(STATION),
DEPTH = mean(DepthBeg),
VSPEED = mean(VSPEED),
TowDistTrack= mean(TowDistTrack),
NetHeight = mean(NetHeight),
NetWidth = mean(NetWidth),
CableOut_S = mean(CableOut_S),
DoorWidth= mean(DoorWidth))## `summarise()` has grouped output by 'Restrictor'. You can override using the
## `.groups` argument.Plotting Aggregate Weights By Species
Looking at some of the more common species. Black line is a 1:1. Models are basic GLM with 95% confidence intervals. Pairs of stations with outlying values are noted. These outliers are retained for the analysis but having some impact and potentially could be removed in some cases (e.g., pair F424 a large school of menhaden was caught in the tow without a restrictor and the vessel moved away for the paired set with the restrictor).
#Could do a lot more like this
library(ggrepel)
Catch %>% as.tibble() %>%
inner_join(.,Station) %>%
mutate(PAIR=str_sub(STATION,-3)) %>%
mutate(SEASON=case_when(month(DATE) < 7 ~ 'S',TRUE~'F')) %>%
mutate(PAIR=paste(SEASON,PAIR,sep='')) %>%
#filter(!PAIR %in% c(479,227,457)) %>%
group_by(VIMSCODE,Restrictor,PAIR) %>% summarise(TotWght=sum(TotWght)) %>%
pivot_wider(names_from = c(Restrictor,VIMSCODE),values_from = TotWght,values_fill = 0) %>%
as_tibble() %>%
pivot_longer(values_to = c('TotWght'),names_to = 'REST_SPP',-PAIR) %>%
mutate(Restrictor=str_sub(REST_SPP,1,1),VIMSCODE=str_sub(REST_SPP,3)) %>%
dplyr::select(-REST_SPP) %>%
pivot_wider(names_from = Restrictor,values_from = TotWght,values_fill=0) %>% #filter(N>0,Y>0) %>%
left_join(.,NM_Species_Code_List_Color_Coordinated %>%
dplyr::select(VIMSCODE=SPCODE,COMMON), by = c('VIMSCODE'='VIMSCODE')) %>%
filter(VIMSCODE %in% c(1028,4,1,171,44,6,15)) %>%
left_join(.,station_pair_info) %>% as_tibble() %>%
#mutate(N=pracma::nthroot(N,4),Y=pracma::nthroot(Y,4)) %>%
#mutate(N=log10(N+1),Y=log10(Y+1)) %>%
left_join(.,station_pair_current %>% dplyr::select(PAIR,CURRENT_DIF)) %>%
left_join(.,pair_zen_data %>% dplyr::select(PAIR,zenith)) %>%
ggplot(aes(x=N,y=Y)) + geom_point(aes(colour=COMMON)) + theme(legend.position = 'None') +
stat_smooth(method = 'glm',level = 0.95,aes(colour=COMMON)) + geom_label_repel(aes(label=PAIR)) +
geom_abline(slope = 1,intercept = 0) + facet_wrap(.~COMMON,scales='free') +
theme(aspect.ratio = 1) + scale_color_tableau()## Joining with `by = join_by(STATION)`
## `summarise()` has grouped output by 'VIMSCODE', 'Restrictor'. You can override
## using the `.groups` argument.
## Joining with `by = join_by(PAIR)`
## Joining with `by = join_by(PAIR)`
## Joining with `by = join_by(PAIR)`
## `geom_smooth()` using formula = 'y ~ x'
Aggregate weight model analysis
This is just one way to explore how well the values from paired stations align. We’re log 10 + constant transforming the data then fitting a simple model to the data. Values on the x axis are catches without the restrictor and values on the y are catches with the restrictor. The plots aren’t exactly what the model is showing, but they are similar.
library(performance)
library(lme4)
library(jtools)
library(sjPlot)
#Making data for fitting an aggregate model
mod_data <- catch_rest %>%
filter(VIMSCODE %in% c(1028,4,1,171,44,6,15)) %>%
left_join(.,station_pair_info) %>%
dplyr::select(-DATE) %>%
mutate(COMMON=as.factor(COMMON),CURRENT_DIF=as.factor(CURRENT_DIF))
#mod_data$COMMON <- fct_rev(mod_data$COMMON)
#mod_data$CURRENT_DIF <- fct_rev(mod_data$CURRENT_DIF)
library(mgcv)
library(gratia)
model_0 <- gam(Y ~ s(N),
data = mod_data, method='REML')
model_1 <- gam(Y ~ s(N)+ s(N,COMMON,bs='sz'),
data = mod_data, method='REML')
model_2 <- gam(Y ~ s(N)+ s(N, by = COMMON),
data = mod_data, method='REML')
compare_performance(model_0,model_1,model_2,rank=TRUE)## # Comparison of Model Performance Indices
##
## Name | Model | R2 | RMSE | Sigma | AIC weights | AICc weights | BIC weights | Performance-Score
## ------------------------------------------------------------------------------------------------------
## model_2 | gam | 0.837 | 0.289 | 0.312 | 0.791 | 0.876 | 1.19e-07 | 63.31%
## model_1 | gam | 0.838 | 0.286 | 0.308 | 0.209 | 0.124 | 2.25e-14 | 56.76%
## model_0 | gam | 0.819 | 0.308 | 0.311 | 2.44e-07 | 5.26e-07 | 1.000 | 23.77%AIC(model_1)## [1] 214.183413(model_1)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ s(N) + s(N, COMMON, bs = "sz")
##
## Estimated degrees of freedom:
## 3.15 16.53 total = 20.68
##
## REML score: 128.602197appraise(model_1)
summary(model_1)##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ s(N) + s(N, COMMON, bs = "sz")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.9192223 0.0203138 45.2512 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(N) 3.14563 3.9099 177.01051 < 2.22e-16 ***
## s(N,COMMON) 16.53388 18.5175 3.57254 1.5241e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.838 Deviance explained = 84.4%
## -REML = 128.6 Scale est. = 0.085202 n = 497#Each species relationship
new_data_1 <-mod_data %>% filter(COMMON == unique(mod_data$COMMON)[1])
new_data_2 <-mod_data %>% filter(COMMON == unique(mod_data$COMMON)[2])
new_data_3 <-mod_data %>% filter(COMMON == unique(mod_data$COMMON)[3])
new_data_4 <-mod_data %>% filter(COMMON == unique(mod_data$COMMON)[4])
new_data_5 <-mod_data %>% filter(COMMON == unique(mod_data$COMMON)[5])
new_data_6 <-mod_data %>% filter(COMMON == unique(mod_data$COMMON)[6])
new_data_7 <-mod_data %>% filter(COMMON == unique(mod_data$COMMON)[7])
new_data <- bind_rows(new_data_1,new_data_2,new_data_3,new_data_4,new_data_5,new_data_6,new_data_7)
m1_pred <- fitted_values(model_1, data = new_data, scale = "response")
p <- ggplot(m1_pred, aes(x = 10^N, y = 10^fitted,
colour = COMMON,group=COMMON)) +
geom_rug(aes(x=10^N,colour=COMMON),sides='b') +
#geom_point(data=aes())
geom_line() +
#geom_point()
geom_ribbon(aes(x=10^N,ymax=10^upper,ymin=10^lower,fill=COMMON,group=COMMON),alpha=0.2) +
facet_wrap(~ COMMON,scales='free',ncol=4) +
geom_abline(slope = 1,intercept = 0) +
labs(color='Species',fill='Species',x='Control Catch (kgs)',
y='Treatment Catch (kgs)') +
scale_color_tableau() + scale_fill_tableau() +
scale_x_log10() + scale_y_log10()
library(lemon)
shift_legend2 <- function(p) {
# ...
# to grob
gp <- ggplotGrob(p)
facet.panels <- grep("^panel", gp[["layout"]][["name"]])
empty.facet.panels <- sapply(facet.panels, function(i) "zeroGrob" %in% class(gp[["grobs"]][[i]]))
empty.facet.panels <- facet.panels[empty.facet.panels]
# establish name of empty panels
empty.facet.panels <- gp[["layout"]][empty.facet.panels, ]
names <- empty.facet.panels$name
# example of names:
#[1] "panel-3-2" "panel-3-3"
# now we just need a simple call to reposition the legend
reposition_legend(p, 'center', panel=names)
}
h <- shift_legend2(p)
#All species combined
new_data_all <- tidyr::expand(mod_data,
N = unique(N))
m1_pred_all <- fitted_values(model_1, data = new_data, scale = "response",exclude='s(N,COMMON)')
g <- ggplot(m1_pred_all, aes(x = 10^N, y = 10^fitted)) +
geom_rug(aes(x=10^N),sides='b') +
#geom_point(data=aes())
geom_line() +
#geom_point()
geom_ribbon(aes(x=10^N,ymax=10^upper,ymin=10^lower),alpha=0.2) +
geom_abline(slope = 1,intercept = 0) +
labs(x='Control Catch (kgs)',
y='Treatment Catch (kgs)') +
scale_color_tableau() + scale_fill_tableau() +
scale_x_log10() + scale_y_log10()
aggregate_weight <- (g + h)
# Making a table of some of the weight model stats
model_fam <- c(model_0 %>% summary() %>% .$family %>% .[1],
model_1 %>% summary() %>% .$family %>% .[1],
model_2 %>% summary() %>% .$family %>% .[1])
model_link <- c(model_0 %>% summary() %>% .$family %>% .[2],
model_1 %>% summary() %>% .$family %>% .[2],
model_2 %>% summary() %>% .$family %>% .[2])
model_formula <- c(model_0 %>% summary() %>% .$formula,
model_1 %>% summary() %>% .$formula,
model_2 %>% summary() %>% .$formula)
model_tab <- tibble(
'Model' = c('Model 0','Model 1','Model 2'),
'Family' = model_fam,
'Link' = model_link,
'Formula' = model_formula,
'n' = c(summary(model_0)$n,summary(model_1)$n,summary(model_2)$n),
'Degrees of freedom' = c(AIC(model_0,model_1,model_2)$df),
'Resid. degrees of freedom' = c(summary(model_0)$residual.df,summary(model_1)$residual.df,summary(model_2)$residual.df),
'R^2' = c(summary(model_0)$r.sq,summary(model_1)$r.sq,summary(model_2)$r.sq),
'Deviance explained' = c(summary(model_0)$dev.expl,summary(model_1)$dev.expl,summary(model_2)$dev.expl),
'AIC' = c(AIC(model_0,model_1,model_2)$AIC)) %>% arrange(`AIC`)Exploring the individual data
Looking at density plots for each of the focus species. The expansion of counts happens using the ‘uncount’ function in dplyr. Definitely curious if that makes sense!
#Making long data for plotting with ggplot
individual_data <- Individual %>%
left_join(.,Station) %>%
filter(VIMSCODE %in% c(1028,4,1,171,44,6,15)) %>%
#uncount(round(ExpFactMeas_SA,0)) %>%
uncount(round(ExpFactMeas,0)) %>%
mutate(SEASON=case_when(month(DATE) < 7 ~ 'Spring',TRUE~'Fall')) %>%
#dplyr::group_by(STATION,Restrictor,VIMSCODE,LENGTH,SEASON) %>%
#dplyr::summarise(count=sum(LENGTH)) %>%
left_join(.,NM_Species_Code_List_Color_Coordinated %>% mutate(VIMSCODE=as.numeric(SPCODE)) %>%
dplyr::select(VIMSCODE,COMMON), by = c('VIMSCODE'='VIMSCODE')) %>%
ungroup() %>%
#removing some outliers
mutate(LENGTH_CM = round(LENGTH/10,0)) %>%
group_by(COMMON) %>%
mutate(LENGTH_NTILE = ntile(LENGTH_CM,n=500)) %>%
filter(LENGTH_NTILE > 1,LENGTH_NTILE < 499) %>%
ungroup()
#Making a plot with mirrored lengths
ml <- ggplot(individual_data %>% mutate(Restrictor = case_when(Restrictor == 'Y' ~ 'Restricted',TRUE ~ 'Standard')), aes(x=LENGTH/10,fill=Restrictor,color=Restrictor,group=Restrictor) ) +
geom_histogram(data = . %>% filter(Restrictor=='Restricted'), aes(y = after_stat(count),fill=Restrictor),binwidth = .25) +
geom_histogram(data = . %>% filter(Restrictor!='Restricted'), aes(y = - after_stat(count),fill=Restrictor),binwidth = .25) +
facet_wrap(~COMMON,scales = 'free',ncol=7) + scale_fill_tableau() + scale_color_tableau() +
labs(color='Tow type',fill="Tow type",x='Length (cm)',y='No individuals')
ml -> mirrored_lengths
mirrored_lengths
Length based modeling
Making data and testing for over/under dispersion
Fitting a regular GLM to the individual data to determine if a binomial model is appropriate (i.e., ‘Dispersion parameter for binomial family taken to be 1’)
#Making wide format data for fitting binomial model
length_mod_data <- Individual %>% filter(VIMSCODE %in% c(1028,4,1,171,44,6,15)) %>%
#::skim()
left_join(.,Station) %>%
#nrow()
uncount(round(ExpFactMeas,0)) %>%
mutate(PAIR=str_sub(STATION,-3)) %>%
mutate(SEASON=case_when(month(DATE) < 7 ~ 'S',TRUE~'F')) %>%
mutate(PAIR=paste(SEASON,PAIR,sep='')) %>%
dplyr::group_by(STATION,PAIR,Restrictor,VIMSCODE,LENGTH,SEASON) %>%
tally() %>% mutate(count = n) %>% dplyr::select(-n) %>%
pivot_wider(values_from = count,names_from = Restrictor,-STATION,values_fill=0) %>%
left_join(.,NM_Species_Code_List_Color_Coordinated %>% mutate(VIMSCODE=as.numeric(SPCODE)) %>%
dplyr::select(VIMSCODE,COMMON), by = c('VIMSCODE'='VIMSCODE')) %>%
left_join(.,station_pair_current %>% dplyr::select(PAIR,CURRENT_DIF)) %>%
left_join(.,pair_zen_data %>% dplyr::select(PAIR,zenith)) %>%
mutate(LENGTH_CM = round(LENGTH/10,0)) %>%
group_by(PAIR,VIMSCODE,LENGTH_CM,SEASON,COMMON,CURRENT_DIF,zenith) %>%
summarise(N=sum(N),Y=sum(Y),TOT = sum(N,Y),PROP=Y/TOT) %>%
group_by(COMMON) %>%
mutate(LENGTH_NTILE = ntile(LENGTH_CM,100)) %>%
filter(LENGTH_NTILE > 10,LENGTH_NTILE < 90)
# Plot with a loess smooth
t <- length_mod_data %>% ggplot(aes(x=LENGTH_CM,y=PROP,colour=COMMON,fill=COMMON)) + geom_point(alpha=0.1) + stat_smooth(method="loess") + facet_wrap(~COMMON,ncol=4,scales='free') + scale_color_tableau() + scale_fill_tableau() + labs(x='Length (cm)',y='Proportion of catch in treatment',colour='Species',fill='Species')
catch_at_length_exploration <- shift_legend2(t)
#Adding averaged station covariates
station_pair_info_avg <- Station %>% mutate(PAIR=str_sub(STATION,-3)) %>%
mutate(SEASON=case_when(month(DATE) < 7 ~ 'S',TRUE~'F')) %>%
mutate(PAIR=paste(SEASON,PAIR,sep='')) %>%
group_by(PAIR) %>%
dplyr::summarise(DoorWidth=mean(DoorWidth),SCOPE=mean(SCOPE),CableOut_P=mean(CableOut_P),
LATITUDE_END=mean(LATITUDE_END),
LONGITUDE_END=mean(LONGITUDE_END),
LATITUDE_START=mean(LATITUDE_START),
LONGITUDE_START=mean(LONGITUDE_START),
NetWidth=mean(NetWidth),
NetHeight=mean(NetHeight),TOWDIST=mean(TOWDIST),
VSPEED=mean(VSPEED),DepthEnd=mean(DepthEnd),DepthBeg=mean(DepthBeg))
#Looking at which tow was first
station_pair_info_order <- Station %>% mutate(PAIR=str_sub(STATION,-3)) %>%
mutate(SEASON=case_when(month(DATE) < 7 ~ 'S',TRUE~'F')) %>%
mutate(PAIR=paste(SEASON,PAIR,sep='')) %>%
dplyr::select(Restrictor,PAIR,TowBegin) %>%
mutate(TowBegin=as.POSIXct(TowBegin,format='%Y-%M-%D %H:%M:%S')) %>%
pivot_wider(names_from = Restrictor,values_from = TowBegin) %>%
mutate(ORDER=case_when(Y > N ~ 'YN',TRUE ~ 'NY')) %>%
dplyr::select(PAIR,ORDER)
#Adding station information
length_mod_data <- length_mod_data %>% left_join(.,station_pair_info_order) %>% left_join(.,station_pair_info_avg)
#Loading mgcv for fitting
library(mgcv)
#Looking at the dispersion in the data
library(dispmod)
length_mod_bi_0 <- glm(cbind(Y,N) ~ LENGTH_CM,
data=length_mod_data %>% mutate(PAIR=as.factor(PAIR)) %>% ungroup(),
family=binomial(link="logit"))
model.disp.test <- glm.binomial.disp(length_mod_bi_0, verbose = T)##
## Binomial overdispersed logit model fitting...
## Iter. 1 phi: 0.242005081
## Iter. 2 phi: 0.278520996
## Iter. 3 phi: 0.278057238
## Iter. 4 phi: 0.278062678
## Iter. 5 phi: 0.278062614
## Converged after 5 iterations.
## Estimated dispersion parameter: 0.278062615
##
## Call:
## glm(formula = cbind(Y, N) ~ LENGTH_CM, family = binomial(link = "logit"),
## data = length_mod_data %>% mutate(PAIR = as.factor(PAIR)) %>%
## ungroup(), weights = disp.weights)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2902676 -1.0097642 0.0141171 0.9761985 2.2306840
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.11118334 0.04553330 2.44180 0.014614 *
## LENGTH_CM -0.00564475 0.00254826 -2.21514 0.026750 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 5681.565 on 4431 degrees of freedom
## Residual deviance: 5676.656 on 4430 degrees of freedom
## AIC: 7144.564
##
## Number of Fisher Scoring iterations: 3print(paste("Dispersion parameter is:", model.disp.test$dispersion))## [1] "Dispersion parameter is: 0.278062614567398"Fitting length-based GAMM models
Here we’re interested in if there are significant effects of the covariates and how large the effects are
#Fitting basic binomial GAM models then adding in random intercepts for the station pair and a random slope for each station pair
library(mgcv)
library(gratia)
length_mod_bi_0 <- gam(PROP ~ s(LENGTH_CM) + s(LENGTH_CM,COMMON,bs='sz') +
s(PAIR,bs='re'),
weights = TOT,
data=length_mod_data %>%
mutate(PAIR=as.factor(as.character(PAIR)),COMMON=as.factor(COMMON)),
family=binomial(link="logit"),method='REML',select=TRUE)
length_mod_bi_1 <- gam(PROP ~ s(LENGTH_CM) + s(LENGTH_CM,COMMON,bs='sz') +
s(PAIR,bs='re') + s(DepthEnd) + s(zenith) +
(ORDER) + (SEASON) + (CURRENT_DIF),
weights = TOT,
data=length_mod_data %>%
mutate(PAIR=as.factor(as.character(PAIR)),COMMON=as.factor(COMMON)),
family=binomial(link="logit"),method='REML',select=TRUE)
length_mod_bi_2 <- gam(PROP ~ s(LENGTH_CM) + s(LENGTH_CM,COMMON,bs='sz') +
s(PAIR,bs='re'),
weights = TOT,
data=length_mod_data %>%
mutate(PAIR=as.factor(as.character(PAIR)),COMMON=as.factor(COMMON)),
family=quasibinomial(link="logit"),method='REML',select=TRUE)
length_mod_bi_3 <- gam(PROP ~ s(LENGTH_CM) + s(LENGTH_CM,COMMON,bs='sz') +
s(PAIR,bs='re') + s(DepthEnd) + s(zenith) +
(ORDER) + (SEASON) + (CURRENT_DIF),
weights = TOT,
data=length_mod_data %>%
mutate(PAIR=as.factor(as.character(PAIR)),COMMON=as.factor(COMMON)),
family=quasibinomial(link="logit"),method='REML',select=TRUE)Comparing model fits
Making plots for the models fit to individual catch data
#Comparing model fit
compare_performance(length_mod_bi_0,length_mod_bi_1,length_mod_bi_2,length_mod_bi_3,metrics = c('R2','RMSE'),rank=TRUE)## # Comparison of Model Performance Indices
##
## Name | Model | R2 | RMSE | Performance-Score
## -----------------------------------------------------------
## length_mod_bi_2 | gam | 0.539 | 0.347 | 100.00%
## length_mod_bi_3 | gam | 0.539 | 0.347 | 73.90%
## length_mod_bi_1 | gam | 0.539 | 0.350 | 23.30%
## length_mod_bi_0 | gam | 0.539 | 0.350 | 0.00%#Making a table of mode lfit statistics
length_model_fam <- c(length_mod_bi_0 %>% summary() %>% .$family %>% .[1],
length_mod_bi_1 %>% summary() %>% .$family %>% .[1],
length_mod_bi_2 %>% summary() %>% .$family %>% .[1],
length_mod_bi_3 %>% summary() %>% .$family %>% .[1]
)
length_model_link <- c(length_mod_bi_0 %>% summary() %>% .$family %>% .[2],
length_mod_bi_1 %>% summary() %>% .$family %>% .[2],
length_mod_bi_2 %>% summary() %>% .$family %>% .[2],
length_mod_bi_3 %>% summary() %>% .$family %>% .[2]
)
length_model_formula <- c(length_mod_bi_0 %>% summary() %>% .$formula,
length_mod_bi_1 %>% summary() %>% .$formula,
length_mod_bi_2 %>% summary() %>% .$formula,
length_mod_bi_3 %>% summary() %>% .$formula
)
length_model_tab <- tibble(
'Model' = c('Model 0','Model 1','Model 2','Model 3'),
'Family' = length_model_fam,
'Link' = length_model_link,
'Formula' = length_model_formula,
'n' = c(summary(length_mod_bi_0)$n,summary(length_mod_bi_1)$n,summary(length_mod_bi_2)$n,summary(length_mod_bi_3)$n),
'Degrees of freedom' = c(AIC(length_mod_bi_0,length_mod_bi_1,length_mod_bi_2,length_mod_bi_3)$df),
'Resid. degrees of freedom' = c(summary(length_mod_bi_0)$residual.df,summary(length_mod_bi_1)$residual.df,summary(length_mod_bi_2)$residual.df,summary(length_mod_bi_3)$residual.df),
'R^2' = c(summary(length_mod_bi_0)$r.sq,summary(length_mod_bi_1)$r.sq,summary(length_mod_bi_2)$r.sq,summary(length_mod_bi_3)$r.sq),
'Deviance explained' = c(summary(length_mod_bi_0)$dev.expl,summary(length_mod_bi_1)$dev.expl,summary(length_mod_bi_2)$dev.expl,summary(length_mod_bi_3)$dev.expl),
'AIC' = c(AIC(length_mod_bi_0,length_mod_bi_1,length_mod_bi_2,length_mod_bi_3)$AIC)) %>% arrange(`Degrees of freedom`)
#Looking at the best model
summary(length_mod_bi_2) ##
## Family: quasibinomial
## Link function: logit
##
## Formula:
## PROP ~ s(LENGTH_CM) + s(LENGTH_CM, COMMON, bs = "sz") + s(PAIR,
## bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0878819 0.1021954 0.85994 0.38987
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(LENGTH_CM) 4.72403 9 11.7944 < 2.22e-16 ***
## s(LENGTH_CM,COMMON) 7.36424 52 106.0900 < 2.22e-16 ***
## s(PAIR) 61.77780 70 58.3801 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.539 Deviance explained = 51.2%
## -REML = 5720 Scale est. = 52.565 n = 4432gam.check(length_mod_bi_2)
##
## Method: REML Optimizer: outer newton
## full convergence after 15 iterations.
## Gradient range [-0.0215448934,0.00951900357]
## (score 5720.04505 & scale 52.5645315).
## Hessian positive definite, eigenvalue range [0.000241508576,2215.94529].
## Model rank = 141 / 141
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(LENGTH_CM) 9.00 4.72 1.25 1
## s(LENGTH_CM,COMMON) 60.00 7.36 NA NA
## s(PAIR) 71.00 61.78 NA NAsummary(length_mod_bi_2)##
## Family: quasibinomial
## Link function: logit
##
## Formula:
## PROP ~ s(LENGTH_CM) + s(LENGTH_CM, COMMON, bs = "sz") + s(PAIR,
## bs = "re")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0878819 0.1021954 0.85994 0.38987
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(LENGTH_CM) 4.72403 9 11.7944 < 2.22e-16 ***
## s(LENGTH_CM,COMMON) 7.36424 52 106.0900 < 2.22e-16 ***
## s(PAIR) 61.77780 70 58.3801 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.539 Deviance explained = 51.2%
## -REML = 5720 Scale est. = 52.565 n = 4432draw(length_mod_bi_2, select = 1) 
pred_df <- NULL
pred_df <- data.frame()
pred_df <- with(length_mod_data,
expand.grid(LENGTH_CM = unique(LENGTH_CM),
COMMON = unique(COMMON)
)
)
# Predictions
pred <- predict(length_mod_bi_2, newdata = pred_df, se.fit = TRUE, type = "link",exclude = c('s(PAIR)'), newdata.guaranteed = TRUE)
pred_df <- cbind(pred_df %>% ungroup(), as.data.frame(pred))
pred_df <- pred_df |> mutate(est = fit,se = se.fit) %>%
mutate(lower = est - (2 * se),
upper = est + (2 * se))
pred_df## LENGTH_CM COMMON fit se.fit
## 1 8 scup 0.300455232158 0.0966166878
## 2 9 scup 0.162453761887 0.0981458249
## 3 10 scup 0.022505277013 0.1033059287
## 4 12 scup -0.146963699692 0.1103857648
## 5 14 scup -0.139416174186 0.1152582378
## 6 15 scup -0.110646947459 0.1198499870
## 7 16 scup -0.089507620581 0.1281003619
## 8 6 scup 0.453679415549 0.1132025132
## 9 7 scup 0.403633681416 0.1013253257
## 10 11 scup -0.086842199961 0.1078506276
## 11 13 scup -0.158901377647 0.1124232858
## 12 23 scup -0.151392176484 0.2605895027
## 13 26 scup -0.228115385999 0.3590068428
## 14 20 scup -0.101021395600 0.1917537203
## 15 22 scup -0.126795426869 0.2351354077
## 16 24 scup -0.181855752630 0.2898990065
## 17 25 scup -0.210389916986 0.3230334918
## 18 27 scup -0.230170055769 0.3974079025
## 19 5 scup 0.454860691481 0.1347364169
## 20 17 scup -0.082184688279 0.1408235557
## 21 29 scup -0.202078085723 0.4823591432
## 22 28 scup -0.218727852695 0.4384848718
## 23 21 scup -0.110839736486 0.2124081764
## 24 18 scup -0.085528509437 0.1563537099
## 25 19 scup -0.093051578822 0.1732262330
## 26 30 scup -0.190411556841 0.5291668269
## 27 32 scup -0.205156521349 0.6299145748
## 28 31 scup -0.190935388417 0.5789209796
## 29 34 scup -0.258942367797 0.7203915087
## 30 36 scup -0.319287031340 0.7900886713
## 31 33 scup -0.229595171316 0.6782240761
## 32 4 scup 0.429390804735 0.1663016399
## 33 35 scup -0.289276167782 0.7565348979
## 34 8 butterfish -0.096104874590 0.0929408075
## 35 9 butterfish -0.207339372085 0.0927701103
## 36 10 butterfish -0.258356888253 0.0934136099
## 37 12 butterfish -0.119476224669 0.0954949564
## 38 14 butterfish 0.139468228122 0.1016553833
## 39 15 butterfish 0.193785165814 0.1063063877
## 40 16 butterfish 0.169682201068 0.1195568585
## 41 6 butterfish 0.148858211489 0.1002649345
## 42 7 butterfish 0.034572987060 0.0940967704
## 43 11 butterfish -0.225954089436 0.0940902730
## 44 13 butterfish 0.020480918078 0.0984508772
## 45 23 butterfish 0.223548409921 0.5392309833
## 46 26 butterfish 0.773853702942 0.7321951915
## 47 20 butterfish -0.080738025198 0.3311681582
## 48 22 butterfish 0.084671820077 0.4736187458
## 49 24 butterfish 0.382458946451 0.6045956886
## 50 25 butterfish 0.562845611066 0.6698935860
## 51 27 butterfish 1.022023520833 0.7891330930
## 52 5 butterfish 0.232960917903 0.1220981457
## 53 17 butterfish 0.087763744108 0.1486951174
## 54 29 butterfish 1.596707497383 0.8900105529
## 55 28 butterfish 1.302443248457 0.8408175682
## 56 21 butterfish -0.023057265443 0.4045472179
## 57 18 butterfish -0.008693597971 0.1958027911
## 58 19 butterfish -0.074424457816 0.2587368499
## 59 30 butterfish 1.878533229417 0.9424086110
## 60 32 butterfish 2.317141785477 1.0757436816
## 61 31 butterfish 2.123443880639 1.0042080661
## 62 34 butterfish 2.558315874789 1.2138967934
## 63 36 butterfish 2.693291950578 1.3150864603
## 64 33 butterfish 2.458725862868 1.1489164033
## 65 4 butterfish 0.296657087860 0.1616933543
## 66 35 butterfish 2.631470314957 1.2676880612
## 67 8 silver hake 0.224389043438 0.1578420533
## 68 9 silver hake 0.122874188933 0.1527090190
## 69 10 silver hake 0.028995322558 0.1484537218
## 70 12 silver hake -0.036077909810 0.1398114536
## 71 14 silver hake 0.057412514681 0.1334707735
## 72 15 silver hake 0.112995546792 0.1315076989
## 73 16 silver hake 0.149414710351 0.1310807116
## 74 6 silver hake 0.336499368402 0.1711967554
## 75 7 silver hake 0.302063643455 0.1639177728
## 76 11 silver hake -0.028444877713 0.1441868106
## 77 13 silver hake -0.000421161808 0.1361802337
## 78 23 silver hake 0.029790602688 0.1520954298
## 79 26 silver hake -0.108063247124 0.1766839843
## 80 20 silver hake 0.125385441251 0.1412415651
## 81 22 silver hake 0.071942463083 0.1485521879
## 82 24 silver hake -0.020175447925 0.1573450409
## 83 25 silver hake -0.069360522422 0.1654718344
## 84 27 silver hake -0.130882167762 0.1918953958
## 85 5 silver hake 0.328968492207 0.1820739939
## 86 17 silver hake 0.162258193778 0.1327377945
## 87 29 silver hake -0.143150799749 0.2379662737
## 88 28 silver hake -0.139818854263 0.2123068427
## 89 21 silver hake 0.103020784979 0.1449944075
## 90 18 silver hake 0.157240658572 0.1354079635
## 91 19 silver hake 0.143200081789 0.1381044602
## 92 30 silver hake -0.150864212044 0.2686430995
## 93 32 silver hake -0.199684793861 0.3422241611
## 94 31 silver hake -0.169531865100 0.3042118287
## 95 34 silver hake -0.275726437753 0.4045257102
## 96 36 silver hake -0.347982532536 0.4415460307
## 97 33 silver hake -0.236947082997 0.3772090809
## 98 4 silver hake 0.297988411913 0.1988939951
## 99 35 silver hake -0.312731967762 0.4244218459
## 100 8 little skate 0.260686228866 0.1004892318
## 101 9 little skate 0.157356310137 0.1005148252
## 102 10 little skate 0.062066870755 0.1014663069
## 103 12 little skate -0.004907852728 0.1019829059
## 104 14 little skate 0.087309302507 0.1044103144
## 105 15 little skate 0.142356652928 0.1067701443
## 106 16 little skate 0.178305836473 0.1106945799
## 107 6 little skate 0.377617048613 0.1046729465
## 108 7 little skate 0.340591228924 0.1016747053
## 109 11 little skate 0.003554395131 0.1019738699
## 110 13 little skate 0.030069813263 0.1027247505
## 111 23 little skate 0.058284434662 0.1501763677
## 112 26 little skate -0.078011459067 0.1738416434
## 113 20 little skate 0.153430417929 0.1352479175
## 114 22 little skate 0.100187179963 0.1460143269
## 115 24 little skate 0.008691259539 0.1554542621
## 116 25 little skate -0.039977252658 0.1632077001
## 117 27 little skate -0.100014755303 0.1884232449
## 118 5 little skate 0.372924435968 0.1134479971
## 119 17 little skate 0.190768851343 0.1166477361
## 120 29 little skate -0.110223200842 0.2332918566
## 121 28 little skate -0.107995764875 0.2082278818
## 122 21 little skate 0.131119309336 0.1411129802
## 123 18 little skate 0.185481531384 0.1232362065
## 124 19 little skate 0.171288387517 0.1293113535
## 125 30 little skate -0.116645014145 0.2633472938
## 126 32 little skate -0.162053949599 0.3355295825
## 127 31 little skate -0.133768450184 0.2982521922
## 128 34 little skate -0.233289680712 0.3958916975
## 129 36 little skate -0.299764924309 0.4302721851
## 130 33 little skate -0.197084969028 0.3696433792
## 131 4 little skate 0.344897931796 0.1312557316
## 132 35 little skate -0.267466696771 0.4145232687
## 133 8 longfin inshore squid 0.615707112072 0.1100071517
## 134 9 longfin inshore squid 0.553064725595 0.1140611131
## 135 10 longfin inshore squid 0.456707249815 0.1204441554
## 136 12 longfin inshore squid 0.284449575047 0.1297356311
## 137 14 longfin inshore squid 0.293352100421 0.1502399417
## 138 15 longfin inshore squid 0.374876674819 0.1630495764
## 139 16 longfin inshore squid 0.485217040762 0.1771827422
## 140 6 longfin inshore squid 0.573357662429 0.1045369440
## 141 7 longfin inshore squid 0.623897470372 0.1082488544
## 142 11 longfin inshore squid 0.357023054037 0.1252214040
## 143 13 longfin inshore squid 0.260691784595 0.1382665120
## 144 23 longfin inshore squid 0.132789175285 0.4742095161
## 145 26 longfin inshore squid -0.748068583869 0.6786750988
## 146 20 longfin inshore squid 0.671672188753 0.2864331393
## 147 22 longfin inshore squid 0.372519017350 0.4066914894
## 148 24 longfin inshore squid -0.144052210575 0.5441660681
## 149 25 longfin inshore squid -0.441376281116 0.6137158888
## 150 27 longfin inshore squid -1.059345728105 0.7367194799
## 151 5 longfin inshore squid 0.476972111591 0.1032519094
## 152 17 longfin inshore squid 0.596817164493 0.1938188663
## 153 29 longfin inshore squid -1.690598812952 0.8368517097
## 154 28 longfin inshore squid -1.374154944324 0.7884786059
## 155 21 longfin inshore squid 0.556935335990 0.3429154278
## 156 18 longfin inshore squid 0.680482761663 0.2141247358
## 157 19 longfin inshore squid 0.711003593426 0.2430938517
## 158 30 longfin inshore squid -2.002007676579 0.8869718295
## 159 32 longfin inshore squid -2.558530464428 1.0088146945
## 160 31 longfin inshore squid -2.296150869960 0.9442762429
## 161 34 longfin inshore squid -2.954589710388 1.1267924919
## 162 36 longfin inshore squid -3.222580792793 1.2034515821
## 163 33 longfin inshore squid -2.778432030145 1.0727554780
## 164 4 longfin inshore squid 0.357764250428 0.1160633462
## 165 35 longfin inshore squid -3.096973853390 1.1685513989
## 166 8 winter flounder 0.260950346357 0.0991889654
## 167 9 winter flounder 0.157608126514 0.0993518909
## 168 10 winter flounder 0.062311016747 0.1004373933
## 169 12 winter flounder -0.004668116940 0.1011689865
## 170 14 winter flounder 0.087551530235 0.1037829473
## 171 15 winter flounder 0.142600498126 0.1062269037
## 172 16 winter flounder 0.178551240286 0.1102314154
## 173 6 winter flounder 0.377918707299 0.1031284363
## 174 7 winter flounder 0.340872214685 0.1002424043
## 175 11 winter flounder 0.003794880065 0.1010603889
## 176 13 winter flounder 0.030310473915 0.1020061225
## 177 23 winter flounder 0.058554267004 0.1500725747
## 178 26 winter flounder -0.077722114578 0.1737573329
## 179 20 winter flounder 0.153686350279 0.1350333910
## 180 22 winter flounder 0.100451954791 0.1458820607
## 181 24 winter flounder 0.008966725999 0.1553688612
## 182 25 winter flounder -0.039695343248 0.1631289476
## 183 27 winter flounder -0.099716921388 0.1883232331
## 184 5 winter flounder 0.373249305530 0.1118730748
## 185 17 winter flounder 0.191015984574 0.1162601780
## 186 29 winter flounder -0.109905208568 0.2331459582
## 187 28 winter flounder -0.107688400344 0.2081060027
## 188 21 winter flounder 0.131379448808 0.1409433552
## 189 18 winter flounder 0.185730891179 0.1229147254
## 190 19 winter flounder 0.171540679419 0.1290463746
## 191 30 winter flounder -0.116315026136 0.2631784025
## 192 32 winter flounder -0.161693723348 0.3353231874
## 193 31 winter flounder -0.133424584830 0.2980633331
## 194 34 winter flounder -0.232888187766 0.3956454676
## 195 36 winter flounder -0.299314156209 0.4299717122
## 196 33 winter flounder -0.196705501525 0.3694188455
## 197 4 winter flounder 0.345247161443 0.1297574574
## 198 35 winter flounder -0.267041083261 0.4142510799
## 199 8 red hake 0.261012365539 0.0993511651
## 200 9 red hake 0.157666087765 0.0994988541
## 201 10 red hake 0.062366139068 0.1005688740
## 202 12 red hake -0.004615933089 0.1012755956
## 203 14 red hake 0.087604560469 0.1038674079
## 204 15 red hake 0.142655870692 0.1063008561
## 205 16 red hake 0.178610376100 0.1102949594
## 206 6 red hake 0.377993155799 0.1033132758
## 207 7 red hake 0.340939733492 0.1004179622
## 208 11 red hake 0.003848118589 0.1011784841
## 209 13 red hake 0.030362513689 0.1021015783
## 210 23 red hake 0.058662151594 0.1500818721
## 211 26 red hake -0.077587735054 0.1737511224
## 212 20 red hake 0.153770491929 0.1350625344
## 213 22 red hake 0.100551541950 0.1458975777
## 214 24 red hake 0.009083261759 0.1553723981
## 215 25 red hake -0.039569922678 0.1631273456
## 216 27 red hake -0.099573632090 0.1883125759
## 217 5 red hake 0.373331780859 0.1120559691
## 218 17 red hake 0.191080169224 0.1163135927
## 219 29 red hake -0.109744259812 0.2331257484
## 220 28 red hake -0.107536278987 0.2080907133
## 221 21 red hake 0.131471129700 0.1409654865
## 222 18 red hake 0.185801123210 0.1229591075
## 223 19 red hake 0.171617640397 0.1290828567
## 224 30 red hake -0.116145125805 0.2631533361
## 225 32 red hake -0.161505145867 0.3352914321
## 226 31 red hake -0.133245492413 0.2980342265
## 227 34 red hake -0.232679900467 0.3956124501
## 228 36 red hake -0.299085685256 0.4299413846
## 229 33 red hake -0.196507168405 0.3693858625
## 230 4 red hake 0.345338211172 0.1299256908
## 231 35 red hake -0.266822727209 0.4142190260
## est se lower upper
## 1 0.300455232158 0.0966166878 0.10722185664 0.4936886077
## 2 0.162453761887 0.0981458249 -0.03383788790 0.3587454117
## 3 0.022505277013 0.1033059287 -0.18410658043 0.2291171345
## 4 -0.146963699692 0.1103857648 -0.36773522935 0.0738078300
## 5 -0.139416174186 0.1152582378 -0.36993264970 0.0911003013
## 6 -0.110646947459 0.1198499870 -0.35034692154 0.1290530266
## 7 -0.089507620581 0.1281003619 -0.34570834433 0.1666931032
## 8 0.453679415549 0.1132025132 0.22727438918 0.6800844419
## 9 0.403633681416 0.1013253257 0.20098302996 0.6062843329
## 10 -0.086842199961 0.1078506276 -0.30254345515 0.1288590552
## 11 -0.158901377647 0.1124232858 -0.38374794920 0.0659451939
## 12 -0.151392176484 0.2605895027 -0.67257118189 0.3697868289
## 13 -0.228115385999 0.3590068428 -0.94612907151 0.4898982995
## 14 -0.101021395600 0.1917537203 -0.48452883610 0.2824860449
## 15 -0.126795426869 0.2351354077 -0.59706624230 0.3434753886
## 16 -0.181855752630 0.2898990065 -0.76165376558 0.3979422603
## 17 -0.210389916986 0.3230334918 -0.85645690061 0.4356770666
## 18 -0.230170055769 0.3974079025 -1.02498586073 0.5646457492
## 19 0.454860691481 0.1347364169 0.18538785776 0.7243335252
## 20 -0.082184688279 0.1408235557 -0.36383179968 0.1994624231
## 21 -0.202078085723 0.4823591432 -1.16679637216 0.7626402007
## 22 -0.218727852695 0.4384848718 -1.09569759637 0.6582418910
## 23 -0.110839736486 0.2124081764 -0.53565608927 0.3139766163
## 24 -0.085528509437 0.1563537099 -0.39823592919 0.2271789103
## 25 -0.093051578822 0.1732262330 -0.43950404473 0.2534008871
## 26 -0.190411556841 0.5291668269 -1.24874521058 0.8679220969
## 27 -0.205156521349 0.6299145748 -1.46498567104 1.0546726283
## 28 -0.190935388417 0.5789209796 -1.34877734772 0.9669065709
## 29 -0.258942367797 0.7203915087 -1.69972538524 1.1818406497
## 30 -0.319287031340 0.7900886713 -1.89946437392 1.2608903112
## 31 -0.229595171316 0.6782240761 -1.58604332345 1.1268529808
## 32 0.429390804735 0.1663016399 0.09678752485 0.7619940846
## 33 -0.289276167782 0.7565348979 -1.80234596354 1.2237936280
## 34 -0.096104874590 0.0929408075 -0.28198648968 0.0897767405
## 35 -0.207339372085 0.0927701103 -0.39287959266 -0.0217991515
## 36 -0.258356888253 0.0934136099 -0.44518410810 -0.0715296684
## 37 -0.119476224669 0.0954949564 -0.31046613753 0.0715136882
## 38 0.139468228122 0.1016553833 -0.06384253843 0.3427789947
## 39 0.193785165814 0.1063063877 -0.01882760949 0.4063979411
## 40 0.169682201068 0.1195568585 -0.06943151593 0.4087959181
## 41 0.148858211489 0.1002649345 -0.05167165757 0.3493880805
## 42 0.034572987060 0.0940967704 -0.15362055367 0.2227665278
## 43 -0.225954089436 0.0940902730 -0.41413463535 -0.0377735435
## 44 0.020480918078 0.0984508772 -0.17642083631 0.2173826725
## 45 0.223548409921 0.5392309833 -0.85491355663 1.3020103765
## 46 0.773853702942 0.7321951915 -0.69053667998 2.2382440859
## 47 -0.080738025198 0.3311681582 -0.74307434159 0.5815982912
## 48 0.084671820077 0.4736187458 -0.86256567156 1.0319093117
## 49 0.382458946451 0.6045956886 -0.82673243075 1.5916503237
## 50 0.562845611066 0.6698935860 -0.77694156100 1.9026327831
## 51 1.022023520833 0.7891330930 -0.55624266507 2.6002897067
## 52 0.232960917903 0.1220981457 -0.01123537342 0.4771572092
## 53 0.087763744108 0.1486951174 -0.20962649076 0.3851539790
## 54 1.596707497383 0.8900105529 -0.18331360835 3.3767286031
## 55 1.302443248457 0.8408175682 -0.37919188787 2.9840783848
## 56 -0.023057265443 0.4045472179 -0.83215170121 0.7860371703
## 57 -0.008693597971 0.1958027911 -0.40029918010 0.3829119842
## 58 -0.074424457816 0.2587368499 -0.59189815754 0.4430492419
## 59 1.878533229417 0.9424086110 -0.00628399268 3.7633504515
## 60 2.317141785477 1.0757436816 0.16565442225 4.4686291487
## 61 2.123443880639 1.0042080661 0.11502774838 4.1318600129
## 62 2.558315874789 1.2138967934 0.13052228808 4.9861094615
## 63 2.693291950578 1.3150864603 0.06311902990 5.3234648713
## 64 2.458725862868 1.1489164033 0.16089305633 4.7565586694
## 65 0.296657087860 0.1616933543 -0.02672962065 0.6200437964
## 66 2.631470314957 1.2676880612 0.09609419259 5.1668464373
## 67 0.224389043438 0.1578420533 -0.09129506315 0.5400731500
## 68 0.122874188933 0.1527090190 -0.18254384913 0.4282922270
## 69 0.028995322558 0.1484537218 -0.26791212104 0.3259027662
## 70 -0.036077909810 0.1398114536 -0.31570081707 0.2435449974
## 71 0.057412514681 0.1334707735 -0.20952903241 0.3243540618
## 72 0.112995546792 0.1315076989 -0.15001985106 0.3760109446
## 73 0.149414710351 0.1310807116 -0.11274671288 0.4115761336
## 74 0.336499368402 0.1711967554 -0.00589414241 0.6788928792
## 75 0.302063643455 0.1639177728 -0.02577190213 0.6298991890
## 76 -0.028444877713 0.1441868106 -0.31681849899 0.2599287436
## 77 -0.000421161808 0.1361802337 -0.27278162923 0.2719393056
## 78 0.029790602688 0.1520954298 -0.27440025694 0.3339814623
## 79 -0.108063247124 0.1766839843 -0.46143121570 0.2453047214
## 80 0.125385441251 0.1412415651 -0.15709768887 0.4078685714
## 81 0.071942463083 0.1485521879 -0.22516191271 0.3690468389
## 82 -0.020175447925 0.1573450409 -0.33486552964 0.2945146338
## 83 -0.069360522422 0.1654718344 -0.40030419125 0.2615831464
## 84 -0.130882167762 0.1918953958 -0.51467295943 0.2529086239
## 85 0.328968492207 0.1820739939 -0.03517949555 0.6931164800
## 86 0.162258193778 0.1327377945 -0.10321739518 0.4277337827
## 87 -0.143150799749 0.2379662737 -0.61908334724 0.3327817477
## 88 -0.139818854263 0.2123068427 -0.56443253964 0.2847948311
## 89 0.103020784979 0.1449944075 -0.18696803000 0.3930096000
## 90 0.157240658572 0.1354079635 -0.11357526838 0.4280565855
## 91 0.143200081789 0.1381044602 -0.13300883867 0.4194090022
## 92 -0.150864212044 0.2686430995 -0.68815041113 0.3864219870
## 93 -0.199684793861 0.3422241611 -0.88413311601 0.4847635283
## 94 -0.169531865100 0.3042118287 -0.77795552257 0.4388917924
## 95 -0.275726437753 0.4045257102 -1.08477785821 0.5333249827
## 96 -0.347982532536 0.4415460307 -1.23107459387 0.5351095288
## 97 -0.236947082997 0.3772090809 -0.99136524484 0.5174710788
## 98 0.297988411913 0.1988939951 -0.09979957834 0.6957764022
## 99 -0.312731967762 0.4244218459 -1.16157565958 0.5361117241
## 100 0.260686228866 0.1004892318 0.05970776532 0.4616646924
## 101 0.157356310137 0.1005148252 -0.04367334029 0.3583859606
## 102 0.062066870755 0.1014663069 -0.14086574307 0.2649994846
## 103 -0.004907852728 0.1019829059 -0.20887366462 0.1990579592
## 104 0.087309302507 0.1044103144 -0.12151132637 0.2961299314
## 105 0.142356652928 0.1067701443 -0.07118363559 0.3558969414
## 106 0.178305836473 0.1106945799 -0.04308332326 0.3996949962
## 107 0.377617048613 0.1046729465 0.16827115557 0.5869629417
## 108 0.340591228924 0.1016747053 0.13724181826 0.5439406396
## 109 0.003554395131 0.1019738699 -0.20039334466 0.2075021349
## 110 0.030069813263 0.1027247505 -0.17537968778 0.2355193143
## 111 0.058284434662 0.1501763677 -0.24206830064 0.3586371700
## 112 -0.078011459067 0.1738416434 -0.42569474577 0.2696718276
## 113 0.153430417929 0.1352479175 -0.11706541707 0.4239262529
## 114 0.100187179963 0.1460143269 -0.19184147391 0.3922158338
## 115 0.008691259539 0.1554542621 -0.30221726472 0.3195997838
## 116 -0.039977252658 0.1632077001 -0.36639265280 0.2864381475
## 117 -0.100014755303 0.1884232449 -0.47686124517 0.2768317346
## 118 0.372924435968 0.1134479971 0.14602844183 0.5998204301
## 119 0.190768851343 0.1166477361 -0.04252662081 0.4240643235
## 120 -0.110223200842 0.2332918566 -0.57680691408 0.3563605124
## 121 -0.107995764875 0.2082278818 -0.52445152839 0.3084599986
## 122 0.131119309336 0.1411129802 -0.15110665104 0.4133452697
## 123 0.185481531384 0.1232362065 -0.06099088171 0.4319539445
## 124 0.171288387517 0.1293113535 -0.08733431945 0.4299110945
## 125 -0.116645014145 0.2633472938 -0.64333960177 0.4100495735
## 126 -0.162053949599 0.3355295825 -0.83311311454 0.5090052153
## 127 -0.133768450184 0.2982521922 -0.73027283452 0.4627359342
## 128 -0.233289680712 0.3958916975 -1.02507307578 0.5584937144
## 129 -0.299764924309 0.4302721851 -1.16030929456 0.5607794459
## 130 -0.197084969028 0.3696433792 -0.93637172741 0.5422017894
## 131 0.344897931796 0.1312557316 0.08238646861 0.6074093950
## 132 -0.267466696771 0.4145232687 -1.09651323424 0.5615798407
## 133 0.615707112072 0.1100071517 0.39569280865 0.8357214155
## 134 0.553064725595 0.1140611131 0.32494249940 0.7811869518
## 135 0.456707249815 0.1204441554 0.21581893899 0.6975955606
## 136 0.284449575047 0.1297356311 0.02497831289 0.5439208372
## 137 0.293352100421 0.1502399417 -0.00712778307 0.5938319839
## 138 0.374876674819 0.1630495764 0.04877752194 0.7009758277
## 139 0.485217040762 0.1771827422 0.13085155636 0.8395825252
## 140 0.573357662429 0.1045369440 0.36428377443 0.7824315504
## 141 0.623897470372 0.1082488544 0.40739976159 0.8403951792
## 142 0.357023054037 0.1252214040 0.10658024611 0.6074658620
## 143 0.260691784595 0.1382665120 -0.01584123932 0.5372248085
## 144 0.132789175285 0.4742095161 -0.81562985692 1.0812082075
## 145 -0.748068583869 0.6786750988 -2.10541878138 0.6092816136
## 146 0.671672188753 0.2864331393 0.09880591011 1.2445384674
## 147 0.372519017350 0.4066914894 -0.44086396148 1.1859019962
## 148 -0.144052210575 0.5441660681 -1.23238434670 0.9442799256
## 149 -0.441376281116 0.6137158888 -1.66880805879 0.7860554966
## 150 -1.059345728105 0.7367194799 -2.53278468800 0.4140932318
## 151 0.476972111591 0.1032519094 0.27046829278 0.6834759304
## 152 0.596817164493 0.1938188663 0.20917943185 0.9844548971
## 153 -1.690598812952 0.8368517097 -3.36430223231 -0.0168953936
## 154 -1.374154944324 0.7884786059 -2.95111215610 0.2028022675
## 155 0.556935335990 0.3429154278 -0.12889551970 1.2427661917
## 156 0.680482761663 0.2141247358 0.25223329003 1.1087322333
## 157 0.711003593426 0.2430938517 0.22481588997 1.1971912969
## 158 -2.002007676579 0.8869718295 -3.77595133559 -0.2280640176
## 159 -2.558530464428 1.0088146945 -4.57615985350 -0.5409010754
## 160 -2.296150869960 0.9442762429 -4.18470335578 -0.4075983841
## 161 -2.954589710388 1.1267924919 -5.20817469411 -0.7010047267
## 162 -3.222580792793 1.2034515821 -5.62948395697 -0.8156776286
## 163 -2.778432030145 1.0727554780 -4.92394298609 -0.6329210742
## 164 0.357764250428 0.1160633462 0.12563755794 0.5898909429
## 165 -3.096973853390 1.1685513989 -5.43407665110 -0.7598710557
## 166 0.260950346357 0.0991889654 0.06257241554 0.4593282772
## 167 0.157608126514 0.0993518909 -0.04109565531 0.3563119083
## 168 0.062311016747 0.1004373933 -0.13856376978 0.2631858033
## 169 -0.004668116940 0.1011689865 -0.20700608997 0.1976698561
## 170 0.087551530235 0.1037829473 -0.12001436437 0.2951174248
## 171 0.142600498126 0.1062269037 -0.06985330933 0.3550543056
## 172 0.178551240286 0.1102314154 -0.04191159053 0.3990140711
## 173 0.377918707299 0.1031284363 0.17166183480 0.5841755798
## 174 0.340872214685 0.1002424043 0.14038740599 0.5413570234
## 175 0.003794880065 0.1010603889 -0.19832589776 0.2059156579
## 176 0.030310473915 0.1020061225 -0.17370177118 0.2343227190
## 177 0.058554267004 0.1500725747 -0.24159088233 0.3586994163
## 178 -0.077722114578 0.1737573329 -0.42523678029 0.2697925511
## 179 0.153686350279 0.1350333910 -0.11638043166 0.4237531322
## 180 0.100451954791 0.1458820607 -0.19131216671 0.3922160763
## 181 0.008966725999 0.1553688612 -0.30177099645 0.3197044484
## 182 -0.039695343248 0.1631289476 -0.36595323854 0.2865625520
## 183 -0.099716921388 0.1883232331 -0.47636338762 0.2769295448
## 184 0.373249305530 0.1118730748 0.14950315586 0.5969954552
## 185 0.191015984574 0.1162601780 -0.04150437138 0.4235363405
## 186 -0.109905208568 0.2331459582 -0.57619712491 0.3563867078
## 187 -0.107688400344 0.2081060027 -0.52390040568 0.3085236050
## 188 0.131379448808 0.1409433552 -0.15050726158 0.4132661592
## 189 0.185730891179 0.1229147254 -0.06009855969 0.4315603420
## 190 0.171540679419 0.1290463746 -0.08655206975 0.4296334286
## 191 -0.116315026136 0.2631784025 -0.64267183123 0.4100417790
## 192 -0.161693723348 0.3353231874 -0.83234009822 0.5089526515
## 193 -0.133424584830 0.2980633331 -0.72955125097 0.4627020813
## 194 -0.232888187766 0.3956454676 -1.02417912294 0.5584027474
## 195 -0.299314156209 0.4299717122 -1.15925758062 0.5606292682
## 196 -0.196705501525 0.3694188455 -0.93554319248 0.5421321894
## 197 0.345247161443 0.1297574574 0.08573224665 0.6047620762
## 198 -0.267041083261 0.4142510799 -1.09554324306 0.5614610765
## 199 0.261012365539 0.0993511651 0.06231003524 0.4597146958
## 200 0.157666087765 0.0994988541 -0.04133162042 0.3566637960
## 201 0.062366139068 0.1005688740 -0.13877160890 0.2635038870
## 202 -0.004615933089 0.1012755956 -0.20716712428 0.1979352581
## 203 0.087604560469 0.1038674079 -0.12013025540 0.2953393763
## 204 0.142655870692 0.1063008561 -0.06994584148 0.3552575829
## 205 0.178610376100 0.1102949594 -0.04197954272 0.3992002949
## 206 0.377993155799 0.1033132758 0.17136660417 0.5846197074
## 207 0.340939733492 0.1004179622 0.14010380912 0.5417756579
## 208 0.003848118589 0.1011784841 -0.19850884955 0.2062050867
## 209 0.030362513689 0.1021015783 -0.17384064285 0.2345656702
## 210 0.058662151594 0.1500818721 -0.24150159262 0.3588258958
## 211 -0.077587735054 0.1737511224 -0.42508997989 0.2699145098
## 212 0.153770491929 0.1350625344 -0.11635457688 0.4238955607
## 213 0.100551541950 0.1458975777 -0.19124361343 0.3923466973
## 214 0.009083261759 0.1553723981 -0.30166153447 0.3198280580
## 215 -0.039569922678 0.1631273456 -0.36582461396 0.2866847686
## 216 -0.099573632090 0.1883125759 -0.47619878391 0.2770515197
## 217 0.373331780859 0.1120559691 0.14921984268 0.5974437190
## 218 0.191080169224 0.1163135927 -0.04154701613 0.4237073546
## 219 -0.109744259812 0.2331257484 -0.57599575669 0.3565072371
## 220 -0.107536278987 0.2080907133 -0.52371770565 0.3086451477
## 221 0.131471129700 0.1409654865 -0.15045984325 0.4134021026
## 222 0.185801123210 0.1229591075 -0.06011709179 0.4317193382
## 223 0.171617640397 0.1290828567 -0.08654807309 0.4297833539
## 224 -0.116145125805 0.2631533361 -0.64245179794 0.4101615463
## 225 -0.161505145867 0.3352914321 -0.83208801005 0.5090777183
## 226 -0.133245492413 0.2980342265 -0.72931394548 0.4628229607
## 227 -0.232679900467 0.3956124501 -1.02390480067 0.5585449997
## 228 -0.299085685256 0.4299413846 -1.15896845444 0.5607970839
## 229 -0.196507168405 0.3693858625 -0.93527889340 0.5422645566
## 230 0.345338211172 0.1299256908 0.08548682964 0.6051895927
## 231 -0.266822727209 0.4142190260 -1.09526077924 0.5616153248pred_df_1 <- pred_df %>% filter(COMMON == unique(length_mod_data$COMMON)[1]) %>% filter(LENGTH_CM > min(length_mod_data %>% filter(COMMON == unique(length_mod_data$COMMON)[1]) %>% .$LENGTH_CM),LENGTH_CM < max(length_mod_data %>% filter(COMMON == unique(length_mod_data$COMMON)[1]) %>% .$LENGTH_CM))
pred_df_2 <- pred_df %>% filter(COMMON == unique(length_mod_data$COMMON)[2]) %>% filter(LENGTH_CM > min(length_mod_data %>% filter(COMMON == unique(length_mod_data$COMMON)[2]) %>% .$LENGTH_CM),LENGTH_CM < max(length_mod_data %>% filter(COMMON == unique(length_mod_data$COMMON)[2]) %>% .$LENGTH_CM))
pred_df_3 <- pred_df %>% filter(COMMON == unique(length_mod_data$COMMON)[3]) %>% filter(LENGTH_CM > min(length_mod_data %>% filter(COMMON == unique(length_mod_data$COMMON)[3]) %>% .$LENGTH_CM),LENGTH_CM < max(length_mod_data %>% filter(COMMON == unique(length_mod_data$COMMON)[3]) %>% .$LENGTH_CM))
pred_df_4 <- pred_df %>% filter(COMMON == unique(length_mod_data$COMMON)[4]) %>% filter(LENGTH_CM > min(length_mod_data %>% filter(COMMON == unique(length_mod_data$COMMON)[4]) %>% .$LENGTH_CM),LENGTH_CM < max(length_mod_data %>% filter(COMMON == unique(length_mod_data$COMMON)[4]) %>% .$LENGTH_CM))
pred_df_5 <- pred_df %>% filter(COMMON == unique(length_mod_data$COMMON)[5]) %>% filter(LENGTH_CM > min(length_mod_data %>% filter(COMMON == unique(length_mod_data$COMMON)[5]) %>% .$LENGTH_CM),LENGTH_CM < max(length_mod_data %>% filter(COMMON == unique(length_mod_data$COMMON)[5]) %>% .$LENGTH_CM))
pred_df_6 <- pred_df %>% filter(COMMON == unique(length_mod_data$COMMON)[6]) %>% filter(LENGTH_CM > min(length_mod_data %>% filter(COMMON == unique(length_mod_data$COMMON)[6]) %>% .$LENGTH_CM),LENGTH_CM < max(length_mod_data %>% filter(COMMON == unique(length_mod_data$COMMON)[6]) %>% .$LENGTH_CM))
pred_df_7 <- pred_df %>% filter(COMMON == unique(length_mod_data$COMMON)[7]) %>% filter(LENGTH_CM > min(length_mod_data %>% filter(COMMON == unique(length_mod_data$COMMON)[7]) %>% .$LENGTH_CM),LENGTH_CM < max(length_mod_data %>% filter(COMMON == unique(length_mod_data$COMMON)[7]) %>% .$LENGTH_CM))
pred_df_new <- bind_rows(pred_df_1,pred_df_2,pred_df_3,pred_df_4,pred_df_5,pred_df_6,pred_df_7) %>% as.tibble() %>% mutate(COMMON=as.character(COMMON))
#Predictions for the ind species
library(boot)##
## Attaching package: 'boot'## The following object is masked from 'package:lattice':
##
## melanomalp <- ggplot(pred_df_new, aes(x = LENGTH_CM, y = inv.logit(est),
colour = COMMON)) +
#geom_rug(aes(x=LENGTH_CM,colour=COMMON),sides='b') +
#geom_point(data=aes())
geom_line() +
#stat_smooth() +
#geom_point() +
geom_ribbon(aes(x=LENGTH_CM,ymax=inv.logit(upper),ymin=inv.logit(lower),fill=COMMON,group=COMMON),alpha=0.2) +
facet_wrap(~ COMMON,scales='free_x',ncol=4) +
geom_hline(yintercept = 0.50) +
scale_y_continuous(limits=c(0,1)) +
scale_color_tableau() + scale_fill_tableau() +
labs(x='Length (cm)',y='Proportion of catch in treatment',colour='Species',fill='Species')
shift_legend2(lp)
#Predictions for the main effect
pred_df_all <- NULL
pred_df_all <- data.frame()
pred_df_all <- with(length_mod_data,
expand.grid(LENGTH_CM = unique(LENGTH_CM)#,
#dum = 0
#COMMON = unique(COMMON)#,
#PAIR = length_mod_data$PAIR[50]
#PAIR = unique(PAIR)
)
)
# Predictions
pred_all <- predict(length_mod_bi_2, newdata = pred_df_all, se.fit = TRUE, type = "link",exclude = c('s(PAIR)','s(LENGTH_CM,COMMON)'),newdata.guaranteed = TRUE)
pred_df_all <- cbind(pred_df_all, as.data.frame(pred_all))
pred_df_all <- pred_df_all |>
mutate(fitted_lower = fit - (2 * se.fit),
fitted_upper = fit + (2 * se.fit))
pred_df_all## LENGTH_CM fit se.fit fitted_lower fitted_upper
## 1 8 0.26101363626 0.0983922816 0.0642290730 0.457798200
## 2 9 0.15766911839 0.0986450097 -0.0396209011 0.354959138
## 3 10 0.06237071253 0.0998169723 -0.1372632320 0.262004657
## 4 12 -0.00460859455 0.1006849194 -0.2059784334 0.196761244
## 5 14 0.08761172318 0.1034134624 -0.1192152016 0.294438648
## 6 15 0.14266049453 0.1059078859 -0.0691552772 0.354476266
## 7 16 0.17861054064 0.1099593611 -0.0413081816 0.398529263
## 8 6 0.37798908137 0.1021685683 0.1736519447 0.582326218
## 9 7 0.34093870849 0.0993581943 0.1422223199 0.539655097
## 10 11 0.00385418296 0.1005136278 -0.1971730727 0.204881439
## 11 13 0.03037042344 0.1015810316 -0.1727916398 0.233532487
## 12 23 0.05860526638 0.1499482775 -0.2412912885 0.358501821
## 13 26 -0.07767354611 0.1736344417 -0.4249424296 0.269595337
## 14 20 0.15374078133 0.1348845362 -0.1160282910 0.423509854
## 15 22 0.10050407862 0.1457547302 -0.1910053817 0.392013539
## 16 24 0.00901668323 0.1552451937 -0.3014737042 0.319507071
## 17 25 -0.03964623258 0.1630054234 -0.3656570794 0.286364614
## 18 27 -0.09966853423 0.1882023519 -0.4760732381 0.276736170
## 19 5 0.37332396222 0.1108896241 0.1515447140 0.595103210
## 20 17 0.19107420275 0.1160308350 -0.0409874673 0.423135873
## 21 29 -0.10985612432 0.2330325138 -0.5759211520 0.356208903
## 22 28 -0.10763983529 0.2079886614 -0.5236171580 0.308337487
## 23 21 0.13143271527 0.1408082045 -0.1501836937 0.413049124
## 24 18 0.18578783694 0.1227201758 -0.0596525147 0.431228189
## 25 19 0.17159633513 0.1288780171 -0.0861596990 0.429352369
## 26 30 -0.11626505459 0.2630679910 -0.6424010366 0.409870927
## 27 32 -0.16164040185 0.3352131637 -0.8320667292 0.508785925
## 28 31 -0.13337325290 0.2979543763 -0.7292820055 0.462535500
## 29 34 -0.23282863001 0.3955216358 -1.0238719017 0.558214642
## 30 36 -0.29924616741 0.4298206271 -1.1588874217 0.560395087
## 31 33 -0.19664943722 0.3693040592 -0.9352575557 0.541958681
## 32 4 0.34532626562 0.1288192144 0.0876878369 0.602964694
## 33 35 -0.26697745446 0.4141145374 -1.0952065293 0.561251620library(boot)
lp_all <- ggplot(pred_df_all, aes(x = LENGTH_CM, y = inv.logit(fit))) +
#geom_rug(aes(x=LENGTH_CM,colour=COMMON),sides='b') +
#geom_point(data=aes())
geom_line() +
#stat_smooth() +
#geom_point() +
geom_ribbon(aes(x=LENGTH_CM,ymax=inv.logit(fitted_upper),ymin=inv.logit(fitted_lower)),alpha=0.2) +
#facet_wrap(~ COMMON,scales='free_x',ncol=4) +
geom_hline(yintercept = 0.50) +
scale_y_continuous(limits=c(0,1)) +
scale_color_tableau() + scale_fill_tableau() +
labs(x='Length (cm)',y='Proportion of catch in treatment',colour='Species',fill='Species')
#Putting things together
lp_all + lp
#sensitivity analysis sand box
#Trying out ggemmeans to reduce coding
#And NTILE lengths to make the length values the same across species (as CMs are different)
#Simpler way to get the same answer
#LENGTH_NTILE model
length_mod_bi_2a <- gam(Y/TOT ~ s(LENGTH_NTILE) + s(LENGTH_NTILE,COMMON,bs='sz') +
s(PAIR,bs='re'),
weights = TOT,
data=length_mod_data %>% #filter(!PAIR %in% c('F479','F547')) %>% without pairs that were large outliers in catch weights
#group_by(PAIR,SEASON) %>% sample_frac(0.75) %>% ungroup() %>% looking at subsets of random pairs
mutate(PAIR=as.factor(as.character(PAIR)),COMMON=as.factor(COMMON)),
family=quasibinomial(link="logit"),method='REML',select=TRUE)
#Comparing to prior models
compare_performance(length_mod_bi_0,length_mod_bi_1,length_mod_bi_2,length_mod_bi_2a,length_mod_bi_3,metrics = c('R2','RMSE'),rank=TRUE)## # Comparison of Model Performance Indices
##
## Name | Model | R2 | RMSE | Performance-Score
## ------------------------------------------------------------
## length_mod_bi_2 | gam | 0.539 | 0.347 | 96.40%
## length_mod_bi_2a | gam | 0.539 | 0.346 | 93.31%
## length_mod_bi_3 | gam | 0.539 | 0.347 | 70.87%
## length_mod_bi_1 | gam | 0.539 | 0.350 | 23.17%
## length_mod_bi_0 | gam | 0.539 | 0.350 | 0.00%#ggeffects giving odd predictions (means too high)
# length_mod_bi_2a %>% ggeffects::predict_response(terms = c('LENGTH_NTILE','COMMON'),margin="mean_reference",type = 'random') %>%
# plot() + facet_wrap(~group) + geom_hline(yintercept = 0.50)
library(ggeffects)
length_mod_bi_2a.emmeans <- ggemmeans(length_mod_bi_2a,
terms=c('LENGTH_NTILE','COMMON'), margin = 'marginalmean',type="count")
plot(length_mod_bi_2a.emmeans, show_data=T) +
labs(title="Predictions by ggemmeans()") + scale_color_tableau() + scale_fill_tableau() +
theme(text = element_text(size = 15)) + facet_wrap(~group) + geom_hline(yintercept = 0.50)## Raw data not available.## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
Looking at gear metrics
library(patchwork)
#Revised gear plots
gear_data <- Station %>% mutate(PAIR=str_sub(STATION,-3)) %>%
mutate(SEASON=case_when(month(DATE) < 7 ~ 'S',TRUE~'F')) %>%
mutate(PAIR=paste(SEASON,PAIR,sep='')) %>%
dplyr::select(STATION,PAIR,SEASON,Restrictor,DoorWidth,SCOPE,CableOut_P,CableOut_S,
LATITUDE_END,NetWidth,NetHeight,TOWDIST,VSPEED,DepthEnd,DepthBeg) %>%
mutate(SEASON=case_when(SEASON == 'S' ~ 'SPRING',TRUE~'FALL'))
gear_data$SEASON <- factor(gear_data$SEASON, levels = c("SPRING", "FALL"))
a <- gear_data %>% #filter(DepthEnd<100) %>%
#mutate(Restrictor = case_when(Restrictor == 'Y' ~ 'Yes',TRUE ~ 'No')) %>%
mutate(Restrictor = case_when(Restrictor == 'N' ~ 'Restricted',TRUE ~ 'Standard')) %>%
filter(NetHeight < 6.25) %>% # Removing an outlier
ggplot(aes(Restrictor,NetHeight,colour=Restrictor)) + geom_boxplot() +
geom_hline(yintercept = 4.7,linetype='dashed') +
geom_hline(yintercept = 5.3,linetype='dashed') +
#stat_summary(color='Black',fun.data = mean_cl_normal, geom = "errorbar", fun.args = list(mult = 1)) +
#geom_point(alpha=0.4) +
theme(legend.position = 'None') +
#geom_label_repel(aes(label=STATION)) +
#facet_wrap(~SEASON) +
scale_x_discrete(limits = rev) +
scale_color_tableau(direction = -1) +
labs(y='Net height (m)',x='')
b <- gear_data %>% #filter(DepthEnd<100) %>%
#mutate(Restrictor = case_when(Restrictor == 'Y' ~ 'Yes',TRUE ~ 'No')) %>%
mutate(Restrictor = case_when(Restrictor == 'Y' ~ 'Restricted',TRUE ~ 'Standard')) %>%
filter(NetWidth>12) %>%
ggplot(aes(Restrictor,NetWidth,colour=Restrictor)) + geom_boxplot() +
geom_hline(yintercept = 14.8,linetype='dashed') +
geom_hline(yintercept = 12.4,linetype='dashed') +
#stat_summary(color='Black',fun.data = mean_cl_normal, geom = "errorbar", fun.args = list(mult = 1)) +
#geom_point(alpha=0.4) +
theme(legend.position = 'None') +
#geom_label_repel(aes(label=STATION)) +
#facet_wrap(~SEASON) +
scale_x_discrete(limits = rev) +
scale_color_tableau(direction = -1) +
labs(y='Net width (m)',x='')
c <- gear_data %>% #filter(DepthEnd<100) %>%
#mutate(Restrictor = case_when(Restrictor == 'Y' ~ 'Yes',TRUE ~ 'No')) %>%
mutate(Restrictor = case_when(Restrictor == 'Y' ~ 'Restricted',TRUE ~ 'Standard')) %>%
ggplot(aes(Restrictor,DoorWidth,colour=Restrictor)) + geom_boxplot() +
#geom_hline(yintercept = 28,linetype='dashed') +
#geom_hline(yintercept = 36,linetype='dashed') +
#stat_summary(color='Black',fun.data = mean_cl_normal, geom = "errorbar", fun.args = list(mult = 1)) +
#geom_point(alpha=0.4) +
theme(legend.position = 'None') +
#geom_label_repel(aes(label=STATION)) +
#facet_wrap(~SEASON) +
scale_x_discrete(limits = rev) +
scale_color_tableau(direction = -1) +
labs(y='Door spread (m)',x='')
d <- gear_data %>%
#mutate(Restrictor = case_when(Restrictor == 'Y' ~ 'Yes',TRUE ~ 'No')) %>%
mutate(Restrictor = case_when(Restrictor == 'Y' ~ 'Restricted',TRUE ~ 'Standard')) %>%
mutate(theta=((DoorWidth * 3.28084)/2)/210,
bridle_angle = asin(theta),bridle_angle = bridle_angle * ( 180.0 / pi )) %>%
ggplot(aes(y=bridle_angle,x=Restrictor,color=Restrictor)) + geom_boxplot() +
scale_color_tableau(direction = -1) +
scale_x_discrete(limits = rev) +
theme(legend.position = 'None') +
labs(y=expression('Bridle angle'~({o})),x='')
gear_metrics_plot <- (a + b + c + d)
gear_metrics_plot
Plotting out products from the script
#Figure 1
sampling_map
#Figure 2
gear_metrics_plot + plot_annotation(tag_levels = 'A')
#Figure 3
catch_ordination + plot_annotation(tag_levels = 'A')
#Figure 4
aggregate_weight + plot_layout(widths = c(1, 2)) + plot_annotation(tag_levels = 'A')
#Figure 5
lp_all + shift_legend2(lp) + plot_layout(widths = c(1, 2)) + plot_annotation(tag_levels = 'A')
#Table 1
station_table %>%
mutate(DEPTH=DEPTH*0.3048) %>%
rename('Season'=SEASON,
'Depth (m)'=DEPTH,
'Vessel speed (kn)'=VSPEED,
'Tow distance (m)'=TowDistTrack,
'Net height (m)'=NetHeight,
'Door width (m)' = DoorWidth,
'Net width (m)'= NetWidth,
'Cable out (ftm)' = CableOut_S) %>%
ungroup() %>% gt() %>% fmt_number(columns=c(4:10),decimals = 2)| Restrictor | Season | Stations | Depth (m) | Vessel speed (kn) | Tow distance (m) | Net height (m) | Net width (m) | Cable out (ftm) | Door width (m) |
|---|---|---|---|---|---|---|---|---|---|
| No | Fall | 31 | 37.89 | 3.07 | 1,837.12 | 5.05 | 14.02 | 95.97 | 33.57 |
| No | Spring | 40 | 36.31 | 3.01 | 1,844.90 | 4.83 | 14.06 | 92.25 | 33.04 |
| Yes | Fall | 31 | 38.07 | 3.06 | 1,902.83 | 5.09 | 13.32 | 95.97 | 31.19 |
| Yes | Spring | 40 | 35.46 | 3.10 | 1,882.66 | 4.91 | 13.63 | 92.38 | 31.64 |
Station %>% mutate(DEPTH=DepthBeg*0.3048) %>% .$DEPTH %>% min()## [1] 19.23288#Table 2
ct| Tot. positive stations | Prop. positive stations | Total weight (kgs) | Station total count | Station mean weight (kgs) | SE station weight (kgs) | CV of station weight (kgs) | Station mean count | SE of station count | CV of station count | |
|---|---|---|---|---|---|---|---|---|---|---|
| longfin inshore squid | ||||||||||
| Y | 71 | 0.97 | 2,734 | 96,787.0 | 38.5 | 7.5 | 2.3 | 1,363.2 | 180.7 | 1.6 |
| N | 68 | 0.93 | 1,739 | 72,938.0 | 24.5 | 2.1 | 1.0 | 1,027.3 | 120.9 | 1.4 |
| little skate | ||||||||||
| Y | 64 | 0.97 | 832 | 1,710.0 | 11.7 | 0.9 | 0.9 | 24.1 | 1.9 | 0.9 |
| N | 63 | 0.95 | 910 | 1,849.0 | 12.8 | 1.1 | 1.1 | 26.0 | 2.4 | 1.1 |
| butterfish | ||||||||||
| N | 61 | 0.97 | 8,140 | 721,698.0 | 114.6 | 13.9 | 1.4 | 10,164.8 | 1,844.2 | 2.2 |
| Y | 60 | 0.95 | 11,988 | 1,015,687.0 | 168.8 | 29.8 | 2.1 | 14,305.5 | 3,746.2 | 3.1 |
| scup | ||||||||||
| N | 53 | 0.96 | 3,074 | 214,630.0 | 43.3 | 6.1 | 1.7 | 3,023.0 | 545.1 | 2.1 |
| Y | 48 | 0.87 | 3,178 | 259,765.0 | 44.8 | 5.6 | 1.5 | 3,658.7 | 630.2 | 2.1 |
| winter flounder | ||||||||||
| Y | 47 | 0.96 | 427 | 2,374.0 | 6.0 | 0.8 | 1.6 | 33.4 | 4.4 | 1.6 |
| N | 43 | 0.88 | 413 | 2,303.0 | 5.8 | 0.8 | 1.6 | 32.4 | 3.9 | 1.4 |
| silver hake | ||||||||||
| Y | 38 | 0.95 | 548 | 6,921.0 | 7.7 | 1.9 | 2.9 | 97.5 | 23.7 | 2.9 |
| N | 37 | 0.92 | 600 | 7,288.0 | 8.4 | 1.4 | 2.0 | 102.6 | 16.2 | 1.9 |
| red hake | ||||||||||
| Y | 28 | 0.93 | 214 | 1,609.0 | 3.0 | 0.7 | 2.6 | 22.7 | 5.8 | 3.0 |
#Table 3
model_tab %>% gt() %>% fmt_number(columns=c(6:10),decimals = 2)| Model | Family | Link | Formula | n | Degrees of freedom | Resid. degrees of freedom | R^2 | Deviance explained | AIC |
|---|---|---|---|---|---|---|---|---|---|
| Model 2 | gaussian | identity | ~, Y, s(N) + s(N, by = COMMON) | 497 | 17.70 | 481.79 | 0.84 | 0.84 | 210.50 |
| Model 1 | gaussian | identity | ~, Y, s(N) + s(N, COMMON, bs = "sz") | 497 | 24.43 | 476.32 | 0.84 | 0.84 | 214.18 |
| Model 0 | gaussian | identity | ~, Y, s(N) | 497 | 3.00 | 495.00 | 0.82 | 0.82 | 244.31 |
#Table 4
length_model_tab %>% gt() %>% fmt_number(columns=c(6:10),decimals = 2)| Model | Family | Link | Formula | n | Degrees of freedom | Resid. degrees of freedom | R^2 | Deviance explained | AIC |
|---|---|---|---|---|---|---|---|---|---|
| Model 2 | quasibinomial | logit | ~, PROP, s(LENGTH_CM) + s(LENGTH_CM, COMMON, bs = "sz") + s(PAIR, bs = "re") | 4432 | 81.43 | 4,357.13 | 0.54 | 0.51 | NA |
| Model 3 | quasibinomial | logit | ~, PROP, s(LENGTH_CM) + s(LENGTH_CM, COMMON, bs = "sz") + s(PAIR, bs = "re") + s(DepthEnd) + s(zenith) + (ORDER) + (SEASON) + (CURRENT_DIF) | 4432 | 82.19 | 4,357.40 | 0.54 | 0.51 | NA |
| Model 1 | binomial | logit | ~, PROP, s(LENGTH_CM) + s(LENGTH_CM, COMMON, bs = "sz") + s(PAIR, bs = "re") + s(DepthEnd) + s(zenith) + (ORDER) + (SEASON) + (CURRENT_DIF) | 4432 | 101.89 | 4,332.15 | 0.54 | 0.51 | 293,752.97 |
| Model 0 | binomial | logit | ~, PROP, s(LENGTH_CM) + s(LENGTH_CM, COMMON, bs = "sz") + s(PAIR, bs = "re") | 4432 | 105.47 | 4,329.55 | 0.54 | 0.51 | 293,731.48 |
#Figure S1
raw_weight_plot + ind_sp_raw_weight_plot + plot_annotation(tag_levels = 'A')
#Figure S2
mirrored_lengths
#Figure S3
shift_legend2(t)## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'